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ABSTRACT 


This report presents an analysis of the in-depth response of materials 
exposed to a high temperature environment, and describes a computer program 
based upon the analysis. 

The differential equations for the in-depth response are formulated 
and then cast into a finite difference form implicit in temperature. Three 
pyrolyzing constituents are allowed, with an accurate model of observed py- 
rolysis kinetics . Heat flow in-depth is one-dimensional, but cross-section 
area may vary with depth. 

The program for in-depth response computation may be coupled to a variety 
of boundary conditions. One of the possible boundary conditions is a com- 
plete boundary layer solution described in other reports of this series. 
Another version of the program may be coupled to a general film coefficient 
model of the boundary layer; this boundary condition is described in some de- 
tail in the present report. 

The report concludes with sample solutions generated by the computer 
program. 
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FOREWORD 

The present report is one of a series of six reports, published simul- 
taneously, which describe analyses and computational procedures for: 1) pre- 

diction of the in-depth response of charring ablation materials, based on one- 
dimensional thermal streamtubes of arbitrary cross-section and considering 
general surface chemical and energy balances, and 2) nonsimilar solution of 
chemically reacting laminar boundary layers, with an approximate formulation 
for unequal diffusion and thermal diffusion coefficients for all species and 
with a general approach to the thermochemical solution of mixed equilibrium- 
nonequilibrium, homogeneous or heterogeneous systems. Part I serves as a 
summary report and describes a procedure for coupling the charring ablator 
and boundary layer routines. The charring ablator procedure is described in 
Part IX, whereas the fluid-mechanical aspects of the boundary layer and the 
boundary-layer solution procedure are treated in Part III- The approximations 
for multicomponent transport properties and the chemical state models are 
described in Parts IV and V, respectively. Finally, in Part VI an analysis 
is presented for the in-depth response of charring materials taking into ac- 
count char-density buildup near the surface due to coking reactions in depth. 

The titles in the series are: 

Part I Summary Report: An Analysis of the Coupled Chemically Reacting 

Boundary Layer and Charring Ablator, by R. M. Kendall, E. P. 
Bartlett, R. A. Rindal, and C. B. Moyer. 

Part II Finite Difference Solution for the In-depth Response of Charring 
Materials Considering Surface Chemical and Energy Balances, by 
C. B. Moyer and R. A. Rindal. 

Part III Nonsimilar Solution of the Multicomponent Laminar Boundary Layer 

by an Integral Matrix Method, by E. P. Bartlett and R. M. Kendall. 

Part IV A Unified Approximation for Mixture Transport Properties for Multi- 
component Boundary-Layer Applications, by E. P. Bartlett, R. M. 
Kendall, and R. A. Rindal. 

Part V A General Approach to the Thermochemical Solution of Mixed Equilib- 
rium-Nonequi librium. Homogeneous or Heterogeneous Systems, by 
R. M. Kendall. 

Part VI An Approach for Characterizing Charring Ablator Response with In- 
depth Coking Reactions, by R. A. Rindal. 

This effort was conducted for the Structures and Mechanics Division of 
the Manned Spacecraft Center, National Aeronautics and Space Administration 
under Contract No. NAS9-4599 to Vidya Division of Itek Corporation with Mr. 
Donald M. Curry and Mr. George Strouhal as the NASA Technical Monitors. The 
work was initiated by the present authors while at Vidya and was completed 
by Aerotherm Corporation under subcontract to Vidya (P.0. 8471 V9002) after 
Aerotherm purchased the physical assets of the Vidya Thermodynamics Depart- 
ment. Dr. Robert M. Kendall of Aerotherm was the Program Manager and Prin- 
cipal investigator. 


v 




t ivL^c 


Dim PACG CUNK not filmed. 


TABLE OF CONTENTS 


ABSTRACT 1_i -i 

FOREWORD V 

LIST OF FIGURES x 

LIST OF SYMBOLS xi 

1. GENERAL INTRODUCTION 1 

2 . PROBLEM DESCRIPTION AND HISTORICAL ORIENTATION 1 

2.1 General Remarks 1 

2.2 Problem Description 2 

2.3 Problem History 3 

3. ANALYSIS AND COMPUTATIONAL PROCEDURE FOR CHARRING MATERIAL 

RESPONSE 4 

3.1 Problem Formulation 4 

# 3 . 1 . 1 Introduction 4 

3.1.2 Basic Differential Equations 4 

3.1.3 Boundary Conditions 6 

3.2 Finite Difference Development 7 

3.2.1 Introduction 7 

3.2.2 Differencing Philosophy 8 

3.2.3 Transformation of Differential Equations 9 

3. 2. 3.1 General Remarks 9 

3.2. 3.2 Geometrical Considerations 10 

3. 2. 3. 3 Conservation of Mass in Moving Coordinate System 10 

3.2.4 Difference Forms 18 

3. 2.4.1 Mass Equation 18 

3. 2. 4. 1.1 Nodes other than the first or last 18 

3. 2. 4. 1.2 The surface node 22 

3. 2. 4. 1.3 The last ablating node 22 

3. 2. 4. 2 Energy Equation 23 

3.2.4. 2.1 Nodes other than the first or last 23 

3. 2. 4. 2. 2 The surface node 25 

3. 2. 4. 2. 3 The last ablating node 27 

3. 2.4. 2.4 Back-up nodes 28 

3. 2. 4. 2. 5 Last node 28 

3.3 Solution Structure Preparatory to Coupling to the Surface 

Boundary Condition 28 

3.3.1 Tri-diagonal Formulation of the Finite Difference 

Energy Relations 28 

3.3.2 Solution of Mass Relations and Evaluation of Tri-diagonal 

Matrix Elements 29 


vii 



TABLE OF CONTENTS (Continued) 




3.3.3 Reduction of Tri-diagonal Matrix to Surface Energy 

Relation 30 

3.4 Coupling In-Depth Response to Surface Energy Balance 30 

3.4.1 General Form of Energy Relation 30 

3.4.2 Tabular Formulation of Surface Quantities 31 

3.4.3 Solution Procedure for the Surface Energy Balance 32 

3.4.4 Completing the In-depth Solution 33 

3.5 Solution Without Energy Balance 

3.6 Solution With Radiation Input Only and No Recession 34 

4. SOME NOTES ON PROPERTY VALUES 35 

4.1 Introduction 35 

4.2 Densities 35 

4.3 Kinetic Data 35 

4.4 Specific Heats 36 

4.5 Heat of Formation 36 

4.6 Heat of Pyrolysis 36 

4.7 Pyrolysis Gas Enthalpy 37 

4.8 Thermal Conductivity 38 

4.9 Surface Emissivity e 39 

4.10 Heat of Ablation, Heat of Combustion 39 

4.11 Surface Species 

4.12 Conclusion 39 

5. NOTES ON FILM COEFFICIENT MODEL OF THE BOUNDARY LAYER WITH 

HEAT TRANSFER, MASS TRANSFER, AND CHEMICAL REACTION 40 

5.1 Introduction 40 

5.2 General Position, Importance, and History of Film Coefficient 

Models Applied to Mass Transfer with Chemical Reactions 40 

5.2.1 General Remarks 40 

5.2.2 History of Film Coefficient Model for C^ / C H 42 

5.2.3 History of Extension to Unequal Mass Diffusion 

Coefficients 43 

5.3 Discussion of Film Coefficient Expressions 43 

5.3.1 Mass Transfer 43 

5.3.2 Energy Equations 44 

5.4 Use of Film Coefficient Equations to Calculate Surface 

Energy Balance 48 

5.4.1 Basic Aspects 48 

5.4.2 Input and Correction of Heat Transfer Coefficient 49 

5.5 Conclusions 50 

6. PRESSURE DROP 50 

6.1 Pressure Drop Correlation Equation 50 


viii 



TABLE OF CONTENTS (Concluded) 


6.2 Finite Difference Formulation 52 

7. SOME OPERATIONAL DETAILS OF THE AERO THERM CHARRING MATERIAL 

'ABLATION PROGRAM, VERSION 2 53 

7.1 Introcution 53 

7.2 Program Description 53 

7.2.1 General Remarks 53 

7.2.2 Program Objectives 54 

7.2.3 Program Capabilities 54 

7.2.4 Solution Procedure 55 

7.2.5 Output Information 55 

7.2.6 Operational Details 56 

7. 2. 6.1 Storage Requirements 56 

7. 2.6. 2 Running Time 56 

7.3 Sample Problem Solutions 56 

7.3.1 Some Typical Problems 56 

7.3.2 Additional Examples 56 

8. SUMMARY AND CONCLUSION 

8.1 General Remarks 57 

8.2 Experience with the In-depth Solution Routine (CMA) 58 

8.3 Conclusion 58 

REFERENCES 59 

FIGURES 53 

APPENDIX A - Equations for Coefficients A , B n , C n , and 
in In-Depth Energy Equation Array* 

APPENDIX B - Conduction Solution Check-Out of the Charring 
Material Ablation Program 

SUB-APPENDIX B-l - Error Analysis of Steady State Solution for a 

Constant Properties Semi-Infinite Receding Solid 
with Constant Surface Temperature and Constant 
Surface Recession Rate 

APPENDIX C - Study of Alternative Treatments of the Decomposition 
Kinetics Equation in the Charring Material Ablation 
Program 

SUB-APPENDIX C-l - Notes on Problem Specification 

SUB-APPENDIX C-2 - Estimation of Execution Time 

APPEXDIX D - Transformation of the In-Depth Energy Equation 


ix 



list of figures 


V 


1* Geometrical Configuration and Coordinate System Illustration 63 

2. Finite Difference Representation 54 

3. Experimentally Determined Coefficients for Flow Through 

Porous Media 65 

4. Typical Output Page 66 

5. Predicted Degradation Depth Histories, Nylon Phenolic Reentry 67 

Problem 

6 . Predicted Temperature Histories at Surface and Selected 68 

Thermocouple Locations, Nylon Phenolic Reentry Problem 

7 . Typical Machine Plotted Isotherm Depths 69 


x 



LIST OF SYMBOLS 


A 

cross section area for a node 

(ft 3 ) 


a bw 

area of back wall of last node 

(ft 3 ) 


A 

n 

coefficient in Eq. (58) . See also Appendix A. 

(Btu/ft 3 

°F) 

A* 

n 

value of A after first pass of Gauss 

reduction, see Eq. (62) 

(Btu/ft 3 

°F) 

B 

pre-exponential factor, Eq. (3) 

(sec -1 ) 


B 

n 

coefficient in Eq. (58) . See also Appendix A 

(Btu/ft 3 

°F) 

B* 

n 

value of B after first pass of Gauss 

reduction, see Eq. (62) 

(Btu/ft 3 

°F) 

B' 

defined as B 1 + B 1 
c g 

( ) 


B' 

C 

defined as m /p u C lr 
c' H e e M 

( ) 


B g 

defined as m /p u C w 
g / M e e M 

( ) 


C H 

Stanton number for heat transfer (corrected 
for "blowing" , if necessary) 

( ) 


C H 

H o 

Stanton number for heat transfer not 
corrected for blowing 

( ) 


C M 

Stanton number for mass transfer 

( ) 


C n 

coefficient in Eq. (58) . See also Appendix A 

(Btu/ ft 3 sec) 

C P 

specific heat 

(Btu/lb 0 

F) 

D n 

coefficient in Eq. (58) . See also Appendix A 

(Btu/ft 3 ) 


D* 

n 

value of D after first pass of Gauss 

reduction, see Eq. (62) 

(Btu/ft 3 ) 


D 

constant defined by Eq. (69) 

(f t 3 /sec) 



binary diffusion coefficient 

(ft 8 /sec) 


E/R.E^R 

activation energy for decomposition 

(°R) 


F 

radiation view factor 

( ) 


F i' F j 

empirical factors appearing in Eq. (69) 

( ) 


F 3 ,F 4' F 5 

general denotations for functional relation- 
ship, Eqs. (58), (59) 

( ) 



xi 



LIST OF SYMBOLS (Continued) 



KE 

k 

k i 

L 

Le 




K. 




m . 
i 


general denotation for functional relationship 

recovery enthalpy 

enthalpy 

sensible enthalpy measured above tempera- 
ture T, 

D 

enthalpy of species i at temperature 
enthalpy of gases adjacent to the wall 

enthalpy term defined at Eq. (38) 
total number of identifiable species 
number of nodelets per node 
nodelet index 

diffusion flux of element k away from wall 

total number of elements in system 
mass fraction of species i 

mass fraction of element k (regardless of 
molecular configuration 

kinetic energy contribution to recovery 
enthalpy 

thermal conductivity, also element index 


pre-exponential factor for constituent i, 
analogous to B 

total number of nodes 

Lewis number p^. .0 /k 
r 3 P 


system molecular weight / x j_^j 
system molecular weight, pyrolysis gas 


molecular weight of species i 

reaction order for constituent i, see Eq, (10) 


mass flow rate per unit area from the 
surface 


( — ) 
(Btu/lb) 
(Btu/lb) 

(Btu/lb) 

(Btu/lb) 

(Btu/lb) 

( ) 

( ) 

( ) 

( ) 

(lb/ft s sec) 

( ) 

( ) 

( ) 

(Btu/lb) 

(Btu/ft sec°F) 
or ( ) 

(sec -1 ) 

(---) 

( ) 

(lb/lb mole) 
(lb/lb mole) 
(lb/lb mole) 
(---) 

(lb/ft 2 sec) 


xii 



N 

NBM 

NL 

q 

q cond 
q diff 
q rad in 


q rad out 

q* 


r.f. 


R 

s 

s 


T,t 


T ,T. 
w 1 


x 


X 


x i 


y 


LIST OF SYMBOLS (Continued) 

mass flow rate of char per unit surface area 

mass flow rate of pyrolysis gas through a 
node 

mass flow rate of gas out a unit area of 
surface 

total number of nodes 

denotes first node of back-up material 

denotes last node of ablation material 

heat absorption term for decomposition 

rate of energy conduction into solid 
material at surface 

rate of energy input to solid surface by 
diffusional processes in the boundary layer 

rate of energy input to the surface by 
radiation from the boundary layer or from 
outside the boundary layer 

rate of energy radiated away from surface 

rate of energy removal from surface with 
condensed phases 

boundary layer recovery factor 

universal gas constant 

distance from original location of receding 
surface to current surface location 

rate of change of S (surface recession 
rate) 

temperature 

wall (surface) temperature 

reference temperature for viscosity Eq. (82) 

velocity of gases at edge of boundary layer 
gas velocity (see pv) 

coordinate normal to ablating surface, 
fixed to receding surface 

1 ^c \ 

virgin mass fraction, x - — y — - |1 - — j 

mole fraction of species i 

coordinate normal to ablating surface, 
origin fixed in space relative to back wall 


( lb/ ft 3 sec) 
(lb/sec) 

( lb/ft a sec) 

( ) 

( ) 

( ) 

(Btu/lb) 
(Btu/ ft 3 sec) 

(Btu/ft 3 sec) 

(Btu/ ft 3 sec) 

(Btu/ ft 3 sec) 
(Btu/ft 2 sec) 

( ) 

( ft# 

[ lb-mole 0 R 

(ft) 

( ft/sec) 

(°R) 

(°R) 

<°R) 

( ft/sec) 

( ft/sec) 


( ) 

(ft) 


xiii 



LIST OF SYMBOLS (Continued) 


Z. 


1 


Z* 

l 




a 

a 


ki 


£ 

Y 

r 

A 

A0 

6 

G 


e 


P 


e 


w 


0 





(pv) w 


diffusion driving force 


Eq. (67) 



j 


modified defined by Eq. (66) 

Z* for elements k independent of molecular 
configuration, defined by Eq. (65) 


( ) 


( — ) 
( — ) 


Greek 


coefficient in Eq. (81) 

mass fraction of element k in species i 
absorptivity of the wall 
coefficient in Eq. (81) 

constant equal to 2/3 on empirical basis 

volume fraction of resin in plastic, see 
Eq. (4) 

denotes change 

time step in finite difference solution 
nodal thickness 

error in an equation supposed to equal 
zero, departure from zero 

volume fraction of imagined undecomposed 
material in a sample of partially degraded 
material 

emissivity of the wall 
time 

viscosity of pyrolysis gas 
dimensionless factor defined by Eq. (68) 

density 

mean volume capacity defined by Eq. (28) 


(ft -2 ) 

( — ) 

( — ) 

(ft -1 ) 

( — ) 
ft 3 resin 
ft 3 material 

( ) 

( ) 

(ft) 

(various) 


ft 3 virgin 
ft 3 material 

( ) 


(lb/ft-sec) 

( ) 

(lb/ft 3 ) 
(Btu/ft 3 °R) 

(lb/ft 2 sec) 


xiv 



LIST OF SYMBOLS (Continued) 


P 

P 

P 

P 

P 


c 

u 


e e 


u c„ 
e e H 


e U e C M 


P 

P 

P 


g 

p 

r 


a 


char density 

mass flow at boundary layer edge 
heat transfer convective film coefficient 
mass transfer convective film coefficient 
original density of a pyrolyzing component 

density of pyrolysis gas 

density of virgin plastic 

final density of a pyrolyzing component 

Stefan-Boltzmann constant 


(lb/ft 3 ) 

(lb/ft 3 sec) 

(lb/ft 3 sec) 

(lb/ft 3 sec) 

(lb/ft 3 resin 
or lb/ft 3 
reinforcement) 

(lb/ft 3 ) 

(lb/ft 3 ) 

(lb/ft 3 resin 
or lb/ft 3 
reinforcement) 

(Btu/ft 3 sec°R 4 ) 


Subscripts 


A 

denotes 

one pyrolyzing component of resin 

B 

denotes 

second pyrolyzing component of resin 

BW 

denotes 

back wall 

b 

denotes 

back face of ablating material 

C 

denotes 

reinforcement 

c 

denotes 

char 

d 

denotes 

decomposition 

e 

denotes 

boundary layer outer edge 

g 

denotes 

pyrolysis gas 

H 

see C T7 
H 


i# j 

denotes 

any identifiable species: atom, ion, molecule 

i 

denotes 

index 

pyrolyzing component? otherwise, general 

j 

denotes 

nodelet 

k 

denotes 

element 

L 

denotes 

last node 


XV 



LIST OF SYMBOLS (Concluded) 


i' 


M 

NL 

n 

o 

P 

r 

res 

s 

w 

0 

1,2 


o 

* 


' (prime) 


£ 

# 


denotes last ablating node 
node index 

denotes datum state for enthalpy; see also p Q 

denotes virgin plastic 
see p 

r 

denotes "reservoir" temperature with which back 
wall communicates 

denotes "sensible” or "thermal" component of 
enthalpy 

denotes wall, i.e., heated surface 
see C TT 

H o 

steps in iteration in Eq. (61) , see also 
Superscripts 


on sensible enthalpy, denotes base temperature above 
which sensible enthalpy is measured; on total enthalpy, 
denotes temperature at which enthalpy is measured 

denotes datum temperature for enthalpy 
see A n' B *n' D n' Z i' \ 


refers to all appearances of element 
of molecular configuration 


k regardless 


denotes average over time interval except on h (q.v.) 
denotes "new", at 0 + A 0, except for B‘, , B^ (q.v.) 

Special Symbols 


defined as, equal to by definition 
pounds force 


xvi 



FINITE DIFFERENCE SOLUTION FOR THE IN-DEPTH RESPONSE 
OF CHARRING MATERIALS CONSIDERING SURFACE CHEMICAL 
AND ENERGY BALANCES 

SECTION 1 

GENERAL INTRODUCTION 

The present report describes a particular analysis and an associated 
computer program for predicting the transient thermal response of a charring 
or decomposing material exposed to a high temperature environment. For 
generality, the in-depth computation program may be coupled to associated 
programs which provide a heated surface boundary condition in the form of 
heat transfer rates and chemical erosion rates. 

Section 2 below describes the general problem and offers a historical 
survey of solution attempts. Sections 3 and 4 present the details of the in- 
depth solution, and Section 5 describes one alternative treatment of the heated 
surface boundary condition. Section 6 describes a model for pressure drop in 
the char, and Section 7 gives a short description of the actual computer 
program and provides some examples of its use. 

The general purpose of this report is to collect in one place certain 
descriptive background material and a number of mathematical derivations 
pertinent to the computer programs developed during this study. The chief 
program pertinent to this report is version 2 of the Aerotherm Charring 
Material Ablation Program (the CMA program) . Related programs are the 
Equilibrium Surface Thermochemistry Program (EST) and the Aerotherm Chemical 
Equilibrium Program (ACE). Instructions for operating these programs are not 
included in this report; the relevant User's Manuals are cited in Section 6 
below. 


SECTION 2 

PROBLEM DESCRIPTION AND HISTORICAL ORIENTATION 
2.1 GENERAL REMARKS 

The basic problem is to predict the temperature and density histories of 
a thermally decomposing material exposed to some defined environment which 
supplies heat and which may chemically erode the material surface. The chief 
practical example of such a problem is the calculation of the performance of 
thermal insulation in hyperthermal environments, including the design of re- 
entry vehicle sacrificial heat shields and expendable rocket nozzle materials. 
Other practical problems include the combustion or charring of wood, the baking 
or various plastics, the combustion of solid rocket fuel, and in-depth pyrolysis 
reactions of all kinds. 
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The general prediction problem may conveniently be divided into two 
parts: the construction of a scheme for computing the in-depth behavior, and 
the specification of the heated surface boundary condition. The present paper 
is mainly concerned with the first problem, although the second topic is also 
given extensive discussion. It may be noted in passing that for quasi-steady 
ablation problems (constant wall temperature, steady recession rate, invariant 
temperature profile with respect to the moving surface) , the details of the 
in-depth solution are not necessary for determining the surface temperature 
and the recession rate. The transient problem, on the other hand, does 
require a complete in-depth solution, and hence is a much more elaborate 
problem. 


2.2 PROBLEM DESCRIPTION 


The physical problem may be illustrated as follows: 


gas 



heated surface 


char or residue 

pyrolysis zone 


As the material is heated, the original virgin material (or rather one or 
more components of the original composite virgin material) pyrolyzes and 
yields a pyrolysis gas, which percolates away from the pyrolysis zone, and a 
porous residue, which for most materials of interest is a carbonaceous char, 
possibly reinforced with refractory fibers or cloth. 

Superimposed on this basic problem may be a number of even more complex 
events. The pyrolysis gases percolating through the char may undergo further 
chemical reactions among themselves, and may react with the char, either 
eroding it or depositing additional residue upon it ("coking"). The char 
itself may collapse or fragment from mechanical or thermal stresses, and the 
refractory reinforcements may melt or suffer mechanical damage. Finally, 
various constituents of the residue structure may react chemically with each 
other, changing the nature of the char, and various mechanical forces may 
remove material from the surface. 


2 



Despite these complexities, it is found that the "simple physics" 
described by 


virgin plastic char + gas 

underlies a wide range of problems of technical interest, and for a great 
many materials, such as carbon phenolic, graphite phenolic, and wood, 
constitute all the events of interest. Such events as coking, mechanical 
erosion, melting, and subsurface reactions (other than pyrolysis) are less 
common and generally characterize specific problems. 

Therefore in any effort to compute the in-depth response of pyrolyzing 
materials the first order of business is to characterize the heat conduction 
and the primary pyrolysis reaction, which have useful generality. Particular 
details of special char chemical systems can then be superimposed upon this 
general computational scheme as required. The present effort has been 
mostly devoted to the general conduction-pyrolysis problem. The numerical 
details are described in Section 3. 


2 . 3 PROBLEM HISTORY 

In general terms, the in-depth calculation requires the solution of a 
differential energy transport equation of the form (for one space dimension, 
and neglecting gas flows for the moment) 


a 


k 


at\ 

ax | 


at 


p c p *55 



o 


plus an associated decomposition or charring kinetic relation 


- f(p.t) 


This coupled pair of differential equations in general defies analytical 
treatment and requires, instead, some approximate technique of solution. 
Perhaps the first general attack on this problem was published by Bamford, 
Crank, and Malan in 1946 (Reference 1). This paper presented a temperature- 
implicit finite difference method, later elaborated in Reference 2, which 
became known as the Crank -Nicolson method. The second paper also presented 
two methods suitable for differential analyzer or analog computations. 

Applications of ablating or sacrificial insulations for rocket nozzles 
and for re-entry vehicle heat shields in the 1950's naturally stimulated a 
great deal of development work on charring material predictive techniques. 
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References 3-30 provide a representative sampling of that part of the litera- 
ture which is aimed primarily at the analysis of the in-depth response of 
charring materials. (A more complete bibliography may be found in Reference 
31.) 

Many of these analyses, particularly References 9, 15, 18, 22, 23, 24, 

27, and 28, were associated with the development of computer programs for in- 
depth response computation. All of these programs treat the in-depth problem 
in very much the same way, with only relatively unimportant variations in the 
numerical formulation and the treatment of the pyrolysis kinetics. 

The present computation scheme does not differ to any very great extent 
from the best versions of earlier programs. It does feature a very great 
fidelity to the observed physics of charring materials, a difference scheme 
which rigorously conserves mass and energy, an economical implicit finite dif- 
ference formulation, and the capability of handling general two-dimensional 
geometries with the limitation of one-dimensional energy flow. 


SECTION 3 

ANALYSIS AND COMPUTATIONAL PROCEDURE FOR 
CHARRING MATERIAL RESPONSE 


3.1 PROBLEM FORMULATION 

3.1.1 Introduction 

Analysis of a complete transient charring material ablation problem 
necessarily involves a computation of the internal thermal response of the 
charring material. As discussed in Section 2 above, this computation consti- 
tutes only a transient heat conduction calculation including the effects of 
internal thermal decomposition with pyrolysis gas generation, coupled to an 
appropriate set of boundary conditions. The sections below present the funda- 
mental assumptions and equations invovled in the in-depth solution. 

3.1.2 Basic Differential Equations 

For the basic in-depth solution, it is assumed that thermal conduction 
is one-dimensional ? however, the cross-section area (perpendicular to the 
conduction flux) is allowed to vary with depth in an arbitrary manner. This 
corresponds to a thermal stream tube. Furthermore, it is assumed that any 
pyrolysis gases formed are in thermal equilibrium with the char. For the 
present discussion, it is assumed that the pyrolysis gases do not react chemi- 
cally with the char in any way. Thus coking or further chemical erosion are 
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excluded. These effects will be discussed in Section 3.1.4 below. Finally, 
any pyrolysis gas formed is assumed to pass immediately out through the char, 
that is, it has zero residence time in the char. Cracking or other chemical 
reactions involving only the pyrolysis gases may be simulated with an appro- 
priate gas specific heat. 

The one-dimensional energy differential equation for this problem is 
readily formulated as 


BO (P hA ) y 


a_ 




(m h ) 
g g' 


e 


(i) 


where p is the density, k is the thermal conductivity, A is area, h is 
enthalpy, and m the local gas flow rate. 

The conservation of mass equation is 


5A g 

dy 


e 



( 2 ) 


Evaluation of this expression requires a specification of the decomposition 
rate dp/30) . A great amount of laboratory pyrolysis data, such as that 

presented in Ref. 32, suggest that the decomposition rate may be taken as an 
Arrhenius type of expression 


lt)y= - B 6 


-E/RT 


P - P, 


(3) 


and for even greater generality it lias been found useful and sufficient to 
consider up to three differently decomposing constituents (Ref. 32) 

p = r(p A + p fi ) + (l - r)p c (4) 
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where for example (p + p fi ) might be the density of resin (or analogous binder) 
in the ablating material, would be the density of the reinforcement and 

T the volume fraction of resin in the virgin plastic composite. Each density 
P A' P B' anc ^ Pq * ma Y follow an independent decomposition equation of the type 
of Eq. (3) . 

It is possible to handle the decomposition in other ways than by Eq. (3) . 
A popular simplification is to treat density as a function of temperature 
only. An even more drastic simplification converts the virgin material to 
complete char at one particular "charring temperature." Other techniques 
specify some char thickness as a function of time, or of heating rate. All 
of these simplifications are, of course, open to objection*. Equation (3) is 
not only the most realistic physically, but is usually easy to handle in com- 
putation . 

3.1,3 Boundary Conditions 

Suitable boundary and initial conditions for the set of equations (1) 
through (4) may be readily formulated. The boundary conditions at the front 
and back faces of the ablating material are usually surface energy balances. 

Of these, the front or "active" surface boundary condition is the most com- 
plex. It is handled in slightly different ways depending on which boundary 
layer treatment is being coupled to the in-depth reaponse program. 

Basically, the surface energy balance may be pictured as 



where the indicated control volume is fixed to the receding surface. Energy 
fluxes leaving the control volume include conduction into the material, radia- 
tion away from the surface, energy in any flow of condensed phase material 
such as liquid runoff, and gross blowing at the surface. Energy inputs to the 
control volume include radiation in from the boundary layer and enthalpy fluxes 
due to char and pyrolysis gas mass flow rates. The final input in the sketch 
is denoted q^ff . It includes all diffusive energy fluxes from the gas 
boundary layer. If the in-depth response computation is being coupled to an 
exact boundary layer solution, the term q^ ff will be available directly as 
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a single term (which is, of course, a complex function of many boundary layer 
properties) . If, on the other hand, the in-depth response is being coupled 
to a simplified boundary layer scheme, such as a convective film coefficient 
model, then the term q^ff has a rather complicated appearance. Section 
5 below contains a further discussion of this aspect of the total computation. 

For the present, it suffices to note that computation of the surface 
energy balance requires the following information from the in-depth solution: 

a. the instantaneous pyrolysis gas rate delivered from in-depth to 
the surface, m^ 

b. a relation between the surface tmmperature and the rate of energy 

conduction into the material, q , 

^■cond 

With these two pieces of information the surface energy balance then deter- 
mines the char consumption rate m and the surface temperature T . It 
will be useful to keep in mind that, from this point of view, the purpose of 
the in-depth solution at any instant is to provide information about m^ 

and q , (T ) . In some circumstances, of course, it is of interest merely 
cona w 

to specify the heated surface temperature and recession rate. In this case 
no surface energy balance is required. 

It is usually of interest to have only one ablating surface. The back- 
wall or non-ablating wall boundary condition may be modelled with a film 
coefficient heat transfer equation. 

3.2 FINITE DIFFERENCE DEVELOPMENT 
3.2.1 Introduction 

Section 3.1.2 above sets forth the governing differential equations whose 
solution is required to define the interal response of the charring material. 
As in many other problems, however, the differential equations cannot be 
solved in general, and it is necessary instead to solve finite difference 
equations* which model the differential equations, and the analyst hopes, re- 
tain the same mathematical properties as the original differential equations. 

A number of plausible difference equations can be proposed, and without the 
benefit of actual experience it is generally impossible to select any particu- 
lar differencing scheme as superior to any other. In the past few years a 


*It is possible, of course, to use simpler schemes than finite difference 
equations. Reference 33 is a sample of the several integral analysis 
approaches which have been tried. Such techniques are usually of insuf- 
ficient accuracy to be generally useful, however. 
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few general differencing principles have heen made reasonably clear, however, 
so that the analyst is not completely in the dark. The following section 
offers some background on this topic. 

3.2.2 Differencing Philosophy 

This section sets down the general principles upon which the finite dif- 
ferencing of the governing equations is based. These principles have proved 
sound and useful, particularly for complex problems. 

In common with all difference procedures, the area of interest (here, 
the charring material) is divided into a number of small zones, each consid- 
ered to be homogeneous. All derivatives in the governing differential equa- 
tions are then replaced by some difference expression from zone to zone. 

These zones, called nodes, thus provide the basic conceptual structure upon 
which the differencing procedure is based. 

The following principles of differencing and nodal sizing have been 
followed in the programming work: 

(1) The nodes have a fixed size. This avoids the slight additional com- 
putation complexity of shrinking nodes, and more importantly, makes principle 

(2) easier to satisy, in addition to preserving a useful nodal spacing through- 
out the history of a given problem. 

(2) Since the nodes are fixed in size, not all of them can be retained 
if the surface of the material is receding due to chemical or mechanical 
erosion. From time to time a node must be dropped, and experience shows that 
it is much more preferable to drop nodes from the back (non-ablating) face 

of the material rather than from the front face. See, for example, Ref. 34 
for a discussion of this problem. This means that the nodal network is “tied 
to the receding surface," and that material appears to be flowing through the 
nodes. This involves a transformation of differential equations (1) and (2) 
to a moving coordinate system and somewhat complicates the algebra of the 
difference equations modelled on these differential equations. Disposing 
of nodes from the front surface, however, often leads to undesirable oscil- 
lations . 

(3) The difference forms of derivatives are kept simple and are formed 
so as to provide a direct physical analog of the differential event leading 
to the derivative. This approach may be contrasted to those approaches 
which seek elaborate difference approximations to derivative expressions. 
Experience shows that the scheme advocated here, while sometimes at a minor 
disadvantage in accuracy, greatly simplifies the attainment of a major objec- 
tive: a difference scheme which conserves energy and mass. Many of the more 
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elaborate difference schemes fail to meet these " simple" but crucial conser- 
vation criteria, and hence frequently converge to erroneous or spurious solu- 
tions . 

(4) The difference equation for energy is formulated in such a way 
that it reduces to the difference equation for mass conservation when tem- 
peratures and enthalpies are uniform. Any lack of consistency between the 
energy and mass equations complicates, and may entirely defeat, convergence 
to a meaningful result. 

(5) The difference energy equations are written to be " implicit" in 
temperature. That is, all temperatures appearing are taken to be "new" 
unknown temperatures applicable at the end of the current time step. It is 
well established that implicit procedures are generally more economical than 
explicit procedures, at least for the majority of ablation problems of inter- 
est in the current work. 

(6) In constrast to point (5) , the decomposition relations are written 
as "explicit” in temperature. To implicitize temperature in these highly non- 
linear equations necessarily involves either a time-consuming iteration pro- 
cedure or an elaborate linearization, which is not necessary for most 
materials. 

(7) Since experience has shown that material decomposition rates are 
strongly dependent on temperature, it is highly desirable to perform the mass 
balance operations in a different, tighter network than that used for the 
energy balance equations. For greatest generality and utility, the number 

of these mass balance "nodelets" per energy balance "node" should be freely 
selectable . 

3.2.3 Transformation of Differential Equations 
3. 2. 3.1 General Remarks 

Solution of the in-depth response will be by difference equations. The 
most convenient finite difference equations governing particular physical 
phenomena may often be derived directly from considering a control volume 
of finite (not infinitesimal) extent. Such is the case, for example, when 
characterizing the thermochemical response of charring materials having con- 
stant cross-sectional area (the flat plate case) . F6r the present problem, 
however, where the cross-sectional area may be an arbitrary function of dis- 
tance below the surface, substantial simplifications may be realized if the 
equations are first considered in differential form. Specifically, differ- 
ences in cross-sectional area with respect to space and time appear in the 
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difference equations derived from the finite control volume approach. The 
neglect of these terms is difficult to rationalize from the difference equa- 
tions; however, if the differential equation is employed, these terms may 
be demonstrated to vanish identically. The resulting differential equations 
may then be expanded into finite difference form yielding a set of finite 
difference equations substantially simpler than those derived directly from 
analyzing a finite control volume. The differential equations are considered 
first, and subsequently expanded to finite difference form. 

As discussed in Section 3.2.2 above, it is deemed convenient to base 
the difference formulation on a nodal network fixed to the heated surface. 
Since this surface will be receding, material will appear to flow into and 
out of the nodes. The differential equations presented as Eqs. (1) and (2) 
require transformation to a moving coordinate system to include this aspect 
of the problem and to provide the proper model for differencing. The mass 
equation is treated first, in Section 3. 2. 3. 3. Discussion of the energy 
equation follows. First, however, some observations about the geometry are 
required. 

3. 2. 3. 2 Geometrical Considerations 

The generalized geometry being considered and the coordinate system to 
be employed for the subsequent differential and finite difference equations 
are shown in Figure 1. 

Geometrical considerations are introduced to the subsequent equations in 
the form of specification of the cross-sectional area as a function of dis- 
tance below the initial surface, A = A(y) . In the event consideration of 
certain special geometries is desired, for example, cylindrically symmetric 
(A « r°) , these functional relationships may be simply obtained to yield 
the cross-sectional area as a function of distance below the initial surface. 

It is important to observe that the origin of the y coordinate is 
fixed in space (relative, say, to the back wall) . Thus a control volume 
at "constant y" contains a fixed, identified piece of material. The origin 
of the x-coordinate , on the other hand, is tied to the receding heated sur- 
face . 

3. 2. 3. 3 Conservation of Mass In Moving Coordinate System 

Decomposition of the ablation material in-depth is characterized by an 
irreversible reaction of the following form: 

plastic — ► char + gas 
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Initially the ablation material is considered to be all plastic, and after 
decomposition it is all char. Intermediate states need not be considered 
until later. The mass conservation equation is written neglecting the mass 
of the gas at any point as being small compared to the mass of solid material 
and assuming that the transit time of the gas from the point of decomposition 
to the heated surface is small, within these constraints the mass conserva- 
tion equation may be written 



dtf 


< P A) y = 



(5) 


where the subscripts indicate variables held constant when performing partial 
differentiation, but A = A(y) alone, so the term 


as ) 


( 6 ) 


and the mass conservation equation becomes 


!V) . * in'] 

iy J e 38 ) 


(7) 


given earlier as Eq. (2) . In the above equation, m , represents the mass 
flow rate of gas past a point, and its derivative with respect to distance 
is seen to equal the gas generation rate. The total gas flow rate passing 
a point is obtained by integration of Eq. (7) 


Xb 


A le ) d y 

y 


(8) 


The material density, 


and density change rate resulting from decomposition 


(dp/d0)y, is obtained from considering the material formulation. The virgin 
plastic may consist of up to three decomposable constituents, and the density 
of the composite is given by Eq. (4) 


np* + Pr) + u - n p c 


(4) 
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where (* + r ) is the density of the resin, p is the density of the rein- 
forcement, and F is the volume fraction of resin in the virgin plastic com- 
posite. The division of resin into A and B components is a consequence of 
the experimentally observed two-stage decomposition process of phenolic resin. 
The rate of change of density resulting from decomposition is given by differ- 
entiating equation (4) with respect to time at constant y. 



where decomposition of each constituent is given by a rate equation of the 
Arrhenius form. 




for 


i - A , B , C 


( 10 ) 


It is now necessary to relate these density changes at constant y to 
the density changes at constant x since we plan to work with the x co- 
ordinate system. At any instant in time the density may be expressed purely 
as a function of spatial position and time, p = p(y,6) • Then 


Differentiating with respect 


From Figure 1 the x and y 
recession 


d( -' = 3y) e dy + de) y de 


(ID 


to time at constant x yields 


at A = a£A azA + acA 

de J x dy ) Q de ) x de / 


( 12 ) 


coordinates are related to the amount of surface 


from which 


y = S + x 


(13) 


azA 

de; 


dS 

ae 


(14) 
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where the surface recession rate, S, is written as an absolute derivative since 
S = S( : ) alone. Substituting this expression into the above equation yields 



(15) 


whi^h, with Equations (9) and (10), represents the desired form of the mass 
conservation equation. 


3. 2.3.4 Conservation of Energy in Moving Coordinate System 

The energy equation is written first with respect to a spatially fixed 
coordinate system. For this purpose, the following functional relationships 
are presumed: 

h = h(T,p) 

T = T(y,e) 

f = p(y,®) 

A = A (y ) 

S = S (6) 

therefore 

h « h(y,0) 

The differential equation governing the conservation of energy within the char- 
ring material is obtained by considering the control volume in Figure 1 and 
equating the net energy transfer rate to the rate of energy change, The re- 
sulting equation was cited as Equation (1) 


storage conduction convection 



For the numerical solution it is convenient to consider a coordinate system 
fixed to the receding surface, as discussed above. To transform the above 
differential equation, which is written for a point, y = constant, to an 
equation written for the moving coordinate system, x = constant, the storage 
term in Eq. (1) may be related to its counterpart in the moving coordinate 
system by expanding the energy change employing the chain rule; 
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phA = phA (y ,0 ) 


d (phA) = (phA) e dy + (phA) y d6 

Differentiating partially with respect to time at constant x yields: 

|e ( P hA) x = ( P hA, e be\ + ae ( P hA) y < 16 > 

Introducing Equation (14) and rearranging obtains 

to <P hA > y = te ( P hA) x ■ ® a7 ( P hA) e (17) 

Substition of Equation (19) into Equation (1) with the observation that partial 
differentiation with respect to x or y at constant time is equivalent, re- 
sults in the transformed energy equation 



III 


IV 


+ S ^ (phA) e + ^ 


(m h ) 

g g 


(18) 


The above terms will be considered separately below. 

Term I 

te ( P hA) x = P h If) + A h ( P h) x U9) 

X 


It is convenient to express the enthalpy change rate in terms of tempera- 
ture and density change rates. For that purpose it is necessary to set down 
some specific model for computing enthalpy. It is deemed convenient and 
fairly reasonable to imagine that partially pyrolyzed material may be regarded, 
for the purpose of computing material properties, as a mixture of pure char 
and pure plastic. A convenient quantity for algebraic manipulation is then 
e , the volume fraction of imagined undecomposed material in the control vol- 
ume. For undecomposed material c is 1; for pure char e = 0, and for 

ir ir 

intermediate states of decomposition it may be anywhere in between. Then the 
density may be written 


p « e p + (1 - e ) p 

K P r P P 


( 20 ) 
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The total enthalpy per unit volume may be written as the mass weighted aver 
age of the enthalpy of the parts 


ph = £ o h + (1 - € ) p h 

^ P*p p p c 


( 21 ) 


where 


♦/ 


C dT 
P P 


(22) 


and 


j. 

h = h u + [ C dT 

c c J p„ 


(23) 


Differentiating Equation (21) obtains 

de 


dh 


dh 


be ( P h) = Pp h P de~ + Pp e P de^ + p c 


de 


dh 


Pc^c e p^c 


(24) 


Differentiating Equations (22) and (23), and noting that the char and plastic 
heats of formation are constant, we have 


= c £T 

de p de 

P 


and 



ZL 


-p c de 


(25) 


Differentiation of Equation (20) results in 



(26) 


Substitution of Equations (25) and (26) into (24) yields the desired rela- 
tion between enthalpy change rate, temperature change rate, and decomposition 
rate 


a_ 

d6 


(ph) x = 


/ Pp h P - Pc h c A 

V Pp " p c J 




(27) 


where 


pC 4peC 
H P p P P P r 


+ (1 - 


€ ) p C 
p H c p 


(28) 
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The specific heat, C^, is the mass weighted average specific heat of the char 
and virgin plastic parts, and, as such, it represents the specific heat of 
the material evaluated in the absence of chemical reactions. 

Utilizing the above, and Equation (12), Term I in the differential 
Equation (11) may be written as follows: 


h ( p hA) x 


p h (If) + A 

N / x 


Pp h P - p c h c 
Pp - Pc 



(29) 


Term II 

Term II in Equation (18) will not require any modification. 
Term III 

For Term III we have 


S (phA) b = Sph + SA ^ (ph) 


(30) 


Now A = A(y) alone, but y = x + S, and S = S(S) alone, so we may write 
A - A(x,e) 


dA 


dx ♦ &) 


bx J, 


b0 J 


ae 


(31) 


Differentiating partially with respect to time at constant y obtains 




hz . \ 

+ M A 

b0 

J Xx) 

, Xu) 

+ Xe) 


y t 

1 y 

7 X 

since A = A(y) alone, 

bA/bS) y 

= o. 

Also, since 


(32) 


xe /, ae 


-S 


(33) 


Combining the above results in 




3x Jq \ 

1 Substituting Equation (34) into (30) yields a new expression for Term III 


(34) 


® h ( p hA) « = ph If) + 


b0 J 


dx ( P h) e 


(35) 
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Term IV 


Term IV may be written 



( 36 ) 


Since differentiation with respect to x and y at constant time are equi- 
valent, substitution of the mass conservation Equation (7) into the above 
gives 


2_ 

dx 


(m h ) 

g g 



+ h A 

g 



(37) 


Substitution of Equations (25) , (35) , and (37) into the energy differential 

Equation (18) yields 


where 


pc 


dT 

de 


i a_ 

A dx 


kA 





h fe) x + It ( P h) s ® 


(38) 


- A Pp h p ~ P e h n 

Pp " Pc 


The terms in Equation (38) represent, from left to right, the sensible 
energy accumulation, the net conduction, the chemical energy accumulation, 
net energy convected as a consequence of coordinate motion, net energy con- 
vected by the pyrolysis gases passing through, and the energy convected away 
by pyrolysis gases generated at the point. All terms are evaluated per unit 
volume. Note that the pyrolysis gas flow rate past a point, m , is the total 
flow rate (see Equation (8) ) and must be divided by area to be on the same 
basis as the other terms in the equation. 

The finite difference formulation of the above derived differential 
equations for mass conservation (Equations (9), (10) and (15)) and energy 

conseration (Equation (38)) are presented in the two succeeding sections. 

Before those types are considered, however, it is of passing interest to 
note that Equation (38) which turns out to be a convenient form for machine 
treatment, can be cast into more appealing form. Some tedious but straight- 
forward algebra, the details of which are given in Appendix D. yields 
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in dh 


m m y ii 

+ Spc |2 + -=2 ,-2 

F p Sx A ox 


With the energy equation in this form, each term has a more readily perceiv- 
able physical significance. 

As a final note, it may be observed that equations (9) , (10) , and (15) 

for mass conservation and Equation (38) for energy conservation cannot be 
used as models for the last ablating node, which is a shrinking as opposed 
to a drifting node. It will prove not necessary to have differential equa- 
tions for this node; the difference equations may be obtained directly. 

3.2.4 Difference Forms 

3. 2.4.1 Mass Equation 


3. 2. 4. 1.1 Nodes other than the first or last 


Experience with the finite difference solution of the preceding differen- 
tial equations has shown that, generally, a much finer definition of the 
space derivative is required for an accurate solution of the mass conserva- 
tion equation than for the energy conservation equation. This is the case 
because, for most material-boundary condition combinations of interest, the 
density profile through the material is much steeper than the corresponding 
temperature profile. as a result, the finite difference spatial grid size 
selected to represent the mass conservation solution is much smaller than that 
selected to represent the energy equation solution. Space derivatives for 
the energy equation solution are obtained by considering an array of nodes, 
while space derivatives for the mass conservation equations are based upon 
consideration of a number of nodelets in each node. A schematic representa- 
tion of the spatial grid is shown in Figure 2 where it is seen that each 
node (n) remains a fixed distance below the surface (x n ) and has a constant 
thickness, 6 . The last node in the ablation material (n = NL) is the one 
exception to this in that it will continually shrink as surface recession 
proceeds until it vanishes at which point the next to last node becomes the 
new last node. Special treatment afforded the last ablating node to include 
consideration of the shrinking process is described subsequently, in Section 
3. 2. 4. 1.3. Referring to Figure 2 it is noted that each node is subdivided 


into J nodelets, each designated n^ where j ranges from 1 to J for 
each node. The mass balance equation is satisfied for each nodelet in the 


following development. 
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The finite difference representation of Equation (15) results from con- 
sidering the density change rate of a nodelet. Since the nodes and nodelets 

remain a fixed distance below the surface (x . = const) , the density change 

n# 3 

rate of a nodelet corresponds to the time derivative of density at constant x 

p n.j ~ p n,j _ (p n.j-H ~ P n .j )S ( ^ p n.j j (40) 

n » J x d 

where primed quantities are evaluated at the time 0' = 0 + t\Q , and the sub- 
script d refers to decomposition which is the density change rate at con- 
stant y. The nodelet temperature, T . required for evaluating the decom- 

n» 3 

position term (Equation (10)) is obtained by linear interpolation between 
adjacent nodal centers. As a consequence of the subsequent convenience of 
obtaining the total nodal density by simple averaging of nodelet densities, 
the nodelets within a node are selected to be of equal volume. The following 
sketch illustrates the scheme for obtaining nodelet temperatures for the spe- 
cial case of four nodelets per node (J = 4) . 



Nodelet 


From this sketch it is noted, for example, that 


*n , 2 " Vi + <A n 6 n + Ah.^.J/2 L 

or, more generally, 


T - T 

n - 1 


+ i 6 A 

2 + 8 °n n 


n,j T n-i + (A 5 + A 5 )/2 

yJ n n n-i n-i ' 


T - T 
a n-i__ 


A b A & 

-HZ X -L ZX + (j . 0.5) 


(41) 


(42) 


for 


j 




si 
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and 


n,D 


T + TT 
n (A 


T - T 
n+i __n_ 


^ 5 . + A 6 )/2 

n+i n+i n n ' 


l* 


- 0.5) 


0.5 


A 5 
n n 


(43) 


for 



The nodelet density change rate of constituent i resulting from decom- 
position is obtained from Equation (10) utilizing the above nodelet tempera- 
tures . 




(44) 


It is noted that the decomposition rate depends upon two quantities, p. . 

i , n / j 

and T which may vary during a time interval (A0) . If temperature and 

3 

density variations during a time interval are large enough to significantly 
effect the decomposition rate, either the time step size must be reduced or 
a temperature and density more representative of the average over the time 
interval should be employed in order to obtain a stable solution. Most gen- 
erally, for problems of practical interest, a nodelet undergoing decomposi- 
tion is experiencing a density decrease and a temperature increase. Equation 
(44) suggests that using the density, p. . , at the beginning of the time 
interval will result in too large a decomposition rate while using the tem- 
perature, T at the beginning of the time interval will result in too 

u 9 3 

small a decomposition rate, for most cases of interest. Since instabilities 
are usually associated with too large rather than too small a change rate, 
it is appropriate to consider an implicit treatment of the density while 
treating the temperature in an explicit manner, at least as a first try. 

Such an explicit treatment of the temperature is possible because the decom- 
position energy of organic constituents is small and, as such, the coupling 
between energy and mass conservation equations is weak for nodes in the de- 
composition region. 

An effective "implicit" treatment of the density in the decomposition 
equation may be readily obtained from direct integration of Equation (44) 
holding the temperature fixed over the time interval. The following results: 
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( 45 ) 


for i- 1, and 
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(46) 


for rru = 1 

These implicit relations yield a much more stable solution than is ob- 
tained with Equation (44) treating the density explicitly.* The overall 
density change rate of a nodelet resulting from decomposition is obtained by 
summing the decomposition rates of each constituent via Equation (10) . The 
total nodelet density change rate over a time interval due to decomposition 

and coordinate system motion is obtained from Equation (40) . 

Since the energy equation is solved on the basis of full nodes rather 
than nodelets it is necessary to evaluate the total nodal density change 
rate. For an even number of equally sized nodelets in each node, 6^ ^ • A^ ^ 

constant, that is, such that each nodelet has the same volume, the total nodal 

density change is the arithmetic average of the density changes of all of the 
nodelets : 


J 



*Calculational results utilizing both integrated and explicit density depen- 
dence in the decomposition equation are presented subsequently in Appendix B. 
Use of the integrated form was suggested in Reference 1, but this device 
appears to have been largely overlooked in more recent work. 


21 



where it is noted that is equivalent to P n+1 ± • This equation can 

be simplified to 


P 


n 



(47) 


3. 2. 4. 1.2 The surface node 

In order to have the surface node be at the surface temperature, it is 
convenient to consider the first node as a half node, with half as many node- 
lets as the other nodes. With the exception of the surface nodelet, the den- 
sity evaluation for nodelets of the first node can be performed according to 
Equation (43) , just as for all the nodelets discussed so far. 

For the surface nodelet it is necessary for consistency with whatever 
solution procedures are supplying the surface energy balance information that 
any material leaving the surface is pure char and hence has the pure char den- 
sity. If this were not the case, then the surface boundary condition solu- 
tions, which are based upon the idea of pure char injection into the surface 
control volume depicted on page 6, could not be coupled in a consistent way 
with the in-depth response solution.* 

3. 2.4.1. 3 The last ablating node 

The last ablating node must be considered separately, since the rear 
boundary of this node is stationary with respect to a fixed coordinate system. 
Within this node, therefore, a variation occurs between the moving and fixed 
coordinate system. Hence Equation (14) must take the form 



for nodelet j (1 < j < J) and where it is noted that the relative motion of 
the x and y coordinate systems varies from 0 at the rear face of the last 
node to nearly S at its forward face. This results in a slightly modified 
nodelet mass balance equation, namely. 


*This implies that if pyrolysis kinetics are slow a dilemma can arise in the 
form of mass accumulation in the surface nodelet. This rarely happens. When 
it does, the wrong tool is being applied to the problem. 

This restriction on leaving density is only necessary when a chemistry solution 
is being used to provide the surface boundary condition. For specified sur- 
face recession rate the actual surface nodelet density can be used and no 
consistency problems arise. 
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3,2.4. 2 Energy Equation 

3. 2. 4. 2.1 Nodes other than the first or last 

A' finite difference representation of Equation (38) can be formulated in 
a variety of ways. As with the mass equation, however, every effort will be 
made here to preserve a correspondence with a finite nodal energy balance. For 
example, the total enthalpy change rate given by Equation (27) , 

be ( P h) x = P c p ae) + ae ) < 51 > 

X y X 


represents the nodal enthalpy change resulting from a change in density and 
a change in temperature. Since enthalpy is a function of T and p only in 
the present analysis, the path followed in going from one temperature-density 
state to another is inconsequential. A constant temperature path followed by 
a constant density path, as illustrated, 
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yields the following interpretation of the two terms on the right side of 
Equation (51) 
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where h is evaluated at the initial temperature, and P ' c p is evaluated 
at the final density and initial temperature. To be precise should be 

evaluated at some mean temperature as well as at the terminal density, but 
otherwise the relation is exact since h is constant at constant temperature 
and p is constant along the second segment of the selected path. 

Equation (38) can be written in finite difference form as 
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where the use of the nodelet densities in the evaluation of oh in the final 
term is motivated by the desire to maintain a consistency between this equa- 
tion and the mass balances (see Eq. (40)).* In effect the node is considered 
to be at constant temperature but with nodelet to nodelet density variations. 

Explicit temperature treatment of this equation is obvious. Implicit 
temperature treatments of all orders of complexity are also possible and the 
selection of the appropriate treatment must be based on a variety of factors - 
but primarily economy, accuracy, and stability. Because the mass equation is 
solved explicitly with respect to temperature, it is possible to first solve 
the mass equation to obtain decomposition rates over the time interval and 
then, employing these rates, solve the energy equation implicitly in tempera- 
ture as follows: 


*Furthermore, if nodelet densities in the ph term are replaced by node 
densities, the overall sub-surface energy balance will not be independent 
of the choice of enthalpy datum temperature, as can be shown by a lengthy 
algebraic development. 


24 





(54) 


where the primed quantities are evaluated at the time O' = 0 + A0. The 
above equation, with the assumption of constant dh^/dT, C^, and other coef- 
ficients over the time interval for a given node permits the direct solution 
of the energy equation using a simple tri-diagonal matrix solution, as will 
be described below. Note that the constant specific heat, C pn , is evaluated 
at the final density, p^, and the initial temperature, T n# for an implicit 
solution of the energy equation after an explicit solution of the mass equa- 
tion . 

Equation (54) applies for all ablating nodes except the first and last. 
Also, for the second node becomes & 1 since the first node is lo- 

cated at the heated surface. The treatment of the two boundary nodes will 
be considered in the next sections. 


3. 2. 4. 2. 2 The surface node 

As noted in Section 3. 2. 4. 1.2, the surface node is treated as a "half- 
node" with temperature equal to surface temperature. Two characteristics 
of this system are significant: the extremely important surface temperature 

is immediately available for the surface boundary condition specification, 
and errors in energy content of the "constant temperature" nodes tend to 
cancel at successive nodal centers as indicated schematically on the follow- 
ing page. 
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It will be recalled from the introductory Section 3.1.3 on boundary 
conditions that one of the key purposes of the in-depth response solution 
is to provide a function q cond (T w ) . How this is finally accomplished will 
not become clear until Section 3.4 below, but it is clear enough that the 
quantity 3 cond * which ultimately will be calculated as part of the surface 
energy balance, will play the central role in linking the in-depth solution 
to the surface energy balance. 

Therefore the energy input to the first node will be left simply as 
q cond' w iii replace a term of the form 

T n-1 - T n 

TT 2 , W 2 

kA k , A -I 

n n n-1 n-1 

Thus we have the energy difference equation for the first node as 



(equation continued on following page} 
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(equation continued from preceding page) 
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(55) 


3. 2.4. 2. 3 The last ablating node 

The energy balance for this node must also be considered separately. 
Using the mass balance equation above (Eq. (50)) as a guide* the following 
representation of the energy equation for the last node of the ablating mate- 
rial is obtained: 
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where it is noted that the final term is roughly half of its counterpart in 
Equation (53) due to the reduced average coordinate motion of this node. 

This reduction is realized by utilizing the ph difference across only half 
the node rather than the full node. The implicit representation of this 
finite difference equation is: 
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(57) 


*For an isothermal system with all components having the same enthalpy per 
unit mass, the energy equation should reduce to the mass equation if abso- 
lute consistency is achieved. 
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with the assumption of constant coefficients this equation will also fall 
within the framework of the tri-diagonal matrix solution to he discussed 
below. 

The nodelet density in the " leaving" ph term of Equations (56) and (57) 
parallels the nodelet density usage in the mass equation (50) . The "entering" 
ph term retains the nodal density. If this procedure is not followed, a 
lengthy algebraic development shows that the overall in-depth energy solution 
will not be independent of the choice of enthalpy datum state. 

3.2.4. 2.4 Back-up nodes 

Energy equations for the back-up nodes are the same as Equation (54) and 
(57) without the decomposition and convection terms, since the nodal structure 
only moves in the ablating material. The back-up nodal structure remains 
fixed . 

3. 2.4. 2. 5 Last node 

The last node of all, whether an ablating node or a node back-up material, 
does not of course conduct energy to an adjacent node. Hence the conduction 
term 


is replaced by a term 


T - T. 
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where ^ res is a back wall convective heat transfer coefficient (which can 
include the effects of radiation) and T res is some reservoir sink temperature 
with which the last node communicates thermally. 


3.3 SOLUTION STRUCTURE PREPARATORY TO COUPLING TO THE SURFACE 
BOUNDARY CONDITION 


3.3,1 Tri-diagonal Formulation of the Finite Difference Energy Relations 

For a given node n, except the first or last, the finite difference energy 
relation involves three unknown temperature, T^ l , T^ f and T^ +i . For the 
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while 


last node N, there are only two unknown temperature, i and T^ # 

the first node equation involves only T ' and , in addition to the heat 

flux q j . 

^cond 

If we arrange all the energy relations in order we obtain an array of the 

form 


B T ' + 
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F 3 ^ q cond^ 
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A N T N-i + 


b n t n 


F (T ) 
4 res 


(58) 

The expressions for the coefficients A , B , C and D are readily ap- 

n n n n 

apparent for the finite difference energy equations (50) and (51) . For ref- 
erence, the detailed expressions for the coefficients are assembled in Appen- 
dix A. 


3.3.2 Solution of Mass Relations and Evaluation of Tri-diagonal 
Matrix Elements 

It is now possible to see clearly what needs to be done for each time 
step A 0 of the solution in order to prepare for coupling to the surface 
energy balance. First, using the current values of p n and § and T n # 
the mass relations (40) , (44) , (47) , (49) , and (56) can be solved, yielding 

"new" gas flow rates m^. . 

This information may then be used to compute the coefficients of the 
tri-diagonal energy equation matrix. Once this matrix is set up, the required 
surface energy relation ^ conc j = q con d^ T w^ 0 ^ ta ^ ne< ^ directly, as des- 

cribed in the next section. 
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3.3.3 Reduction of Tri-diagonal Matrix to Surface Energy Relation 


Referring to the array of in-depth energy equations set down symbolically 
as Set ( 58 ) in Section 3.3.1 above, it may be seen that, beginning with the 
last node, the highest-indexed unknown temperature may be eliminated from 
each equation of Set (58) in turn. (This is the standard first step in the 
routine reduction of a tri-diagonal matrix.) Of the resulting simpler set of 
equations (shown as Set (62) below and discussed at that time) only the top- 
most is of immediate interest. It may be arranged as 


cond 


F (T‘) 

6 l 


( 59 ) 


where F & is a simple linear relation. It will be noted that this reduction 
implies that the A, B, C, and D terms involve only known quantities evaluated 
at the beginning of the time step. In particular, the surface recession rate 
S is treated in this explicit manner. This causes little error since the 
energy terms in depth involving S are small compared to other energy terms. 

Since T. = T , Equation (59) is the desired relation between q . 

1 w ^cond 

and T w implied by the in-depth solution. 

It is now necessary to harmonize this in-depth relation with the surface 
energy balance. This will be discussed in the following section. 

3.4 COUPLING IN-DEPTH RESPONSE TO SURFACE ENERGY BALANCE 

3.4.1 General Form of Energy Relation 

If the surface boundary condition involves an energy balance with con- 
vective energy input, the final in-depth relation Equation (55) must now be 
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coupled to the surface energy balance illustrated in the sketch above. The 
surface energy balance may be written 

q diff + q rad + "c h c + ™g h g “ q rad “ (pv) w h w ' q * “ q cond = 0 (60) 

in out 

where 

(pv) = m + m - m* 
v H w g g c 

We note that m and q , — F_ (T ) are delivered by the in-depth solution 
g cond 6 w 

Other dependencies of interest are 


h 

= h „( T J ' 

g 

g w 

h 

= h (Tj , 

c 

c w 

q rad 

= q rad ( V 

out 

out 


For the other terms, we have 


m*, q d -ff, q h , q* = functions of boundary layer edge enthalpy, 

° x 1 ra w pressure, boundary layer aerodynamic solu- 

in tion, conservation of chemical element laws, 

chemical equilibria and/or kinetic relations 
upstream events, 


3.4.2 Tabular Formulation of Surface Quantities 

From the standpoint of the surface energy balance solution the desired 
relationship may be summarized as 


m* , 


q dif f ' 


d rad' 

in 


functions of (m ,T ) 
1 c w 


where all the other aspects of the solution are subsumed in some other compu- 
tational procedure, the details of which are not of immediate interest. This 
other computational procedure might be based upon a simple film coefficient 
model of the boundary layer, as will be described below, or it might be based 
on a very detailed boundary layer solution. 
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In almost all cases, however, it proves expedient to do this solution 
outside the charring material solution routine and to make the results avail- 
able to the surface energy balance operation in the form of tables of q^ff, 

q , . r h , q*, and T as functions of m , , and another variable which 

M rad in w M w c 9S 

is essentially time 0 and includes all time dependent aspects such as pres- 
sure . 


3.4.3 Solution Procedure for the Surface Energy Balance 


The energy balance solution procedure is then fairly obvious. An initial 

guess of the char consumption rate, m^, is obtained in some manner. With 

this m c , and with the supplied by the in-depth solution, the quantities 

q q j . , h , q*, and T are obtained by table look up in the tables 

M di f f M rad in' w # M w J ^ 

provided by the outside surface solution routine. The quantities h^, h 

and q . , can then be formulated using the T so obtained. 

H rad out ^ w 

Then the surface energy balance (60) can be computed. In general, however, 
the sum of the terms will not equal zero but some non-zero quantity € called the 
error. Some appropriate iteration procedure must be devised to select succes- 
sively better estimates of m c which drive the error e to zero. Experience 
shows that Newton's procedure, in which the derivative of the error with re- 
spect to m c is used to compute the next guess for m c as schematically il- 
lustrated below, gives good results. 



n = m' 
c c 

2 1 


(ds/dA,) ~ 


( 61 ) 
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Since the actual e function is almost linear between tabular points (radia- 
tion introduces slight curvature) this scheme converges rapidly to an answer. 
Possible traps due to kinks or elbows in the e function are avoided by limit- 
ing corrections on m^ to one tabular interval at a time. 

3.4.4 Completing the In-depth Solution 

Once the surface energy balance has been satisfied to an acceptable level 
of accuracy, the new surface temperature, , may be substituted in the re- 
duced array of temperatures which began as Set (58) . The highest-indexed 
unknown temperature has been eliminated from each equation in this set in the 
process of deriving Equation (59) . The remaining array looks like 


\ T ; + b j t 2 


D 2 


AT' + B*T ' = D* 

3 3 3 3 3 

AT' + B*T ' = D* 

4 3 4 4 4 


A^T + B*T'_ = F (T ) (6 2) 

N N-i N N 4 ' res' v 

S ince T^ is now known, the first equation of Set (62) yields T^ directly, 
then the second equation yields T^ , and so on until the new temperature set 
is complete. 

As a final step, new values for temperature dependent properties can be 
selected for each node and the entire system is then ready for a new time step, 
beginning with the decomposition event. 

Alternatively, the question arises whether it might be prudent to iterate 
the entire time step calculation to implicitize the quantity S (for which a 
new value may be computed from the new m c ) , as well as to implicitize the 
decomposition kinetics and all temperature dependent properties. Usually, 
however, the quantity S exerts only a weak influence on the in-depth solu- 
tion (which tends to adjust itself to damp out the effects of oscillating S 
values) . Furthermore, all important temperature dependent property effects 
have already been partly implicitized in the in-depth energy solution through 
linearization. Thus the S aspect and the properties aspect do not seem to 
provide much motivation for iteration. 
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The decomposition problem is another matter. It is clear that leaving 
the temperature dependence explicit in the decomposition kinetics is inviting 
trouble, since for most practical problems of interest the temperature rise 
is great and the decomposition rates are very sensitively dependent on tem- 
perature. Indeed, practical experience shows that for violent transients 
the problem solution procedure can be disrupted by instabilities due to this 
explicit temperature dependence. 

It is possible that a linearization of this dependence may be sufficient 
to effect the solutions to most problems without the necessity for iteration. 

In that case the in-depth temperature matrix is triangular instead of tri- 
diagonal, but the basic solution procedure is not much modified. This device 
has not yet been incorporated in the charring ablation computer program des- 
cribed in Section 6 below. 

■i . 5 SOLUTION WITHOUT ENERGY BALANCE 

The surface boundary condition need not, of course, be an energy balance. 
Surface temperature and recession rate might be specified. In that case T* 
is known, and the solution of Section 3.4,4 can be completed at once. The 
quantity ^ con( j i s on ^V °f “cultural" interest now. 

This option is especially useful for parametric studies matching inter- 
nal thermocouple response predictions to the measured thermocouple responses, 
using measured surface temperature and recession data, in order to "back-out" 
thermal conductivity data, determining conductivity as a function of tempera- 
ture and density. Frequently this procedure is required before a surface 
energy balance type of calculation can be done since the required thermal 
conductivity information is not generally available for many materials . 

3.6 SOLUTION WITH RADIATION INPUT ONLY AND NO RECESSION 

To simulate "cool-down" or " soak-out" problems, convective energy input 

can be ignored and surface recession supressed. The surface temperature can 

be determined with a simple version of Equation (60) in which the terms 

a * . ^ a „ . , m , and cr* are zero and for which the terms m h„ and 
M diff A rad m c ^ <3 g 

(pv) w h w cancel. In this case the wall temperature is the independent variable. 
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SECTION 4 

SOME NOTES ON PROPERTY VALUES 


4.1 INTRODUCTION 

Some remarks on the various physical property values required for the 
in-depth solution may help to clarify certain aspects of both the analysis 
of Section 3 above and the associated computer program. Features which are 
especially important or somewhat unusual (when compared to most other compu- 
tation schemes) will be stressed. 

4.2 DENSITIES 

Section 3. 2.3.3 above describes how the decomposition kinetic relations 

are written for three independently pyrolyzing components requiring , 

u i 

the initial (virgin) density for each component i, and p r> , the final 
(char) density for each component. Thus the user must supply p Q . and p r _^ 
for each component, and the resin volume fraction T. The program itself 
computes the overall virgin and char densities from this data. 

The user frequently starts with measured values of virgin and char den- 
sities and partitions these densities up among the decomposing constituents 
according to the decomposition kinetic data he is trying to match. The next 
section gives some discussion of this point. 

4.3 KINETIC DATA 

Equation (4b) of Section 3, 2. 3. 3 indicates that in addition to the p Qi 
and p r ^ data, the kinetic constants k^, nu and e^/R are required for 
each constituent of the decomposing material. The required data, much of 
which may be located in the literature, are obtained from thermogravimetric 
(TGA) experiments. The observed decomposition kinetics may be modelled with 
one, two, or three kinetic regimes as required. Experience shows that more 
than three regimes are not required for materials of current interest. 

Indeed, the great majority of materials is well described by only a single 
kinetic relation. 

As noted in Section 4.2 above, the selection of the decomposing consti- 
tuents influences the partition of the densities among the component densi- 
ties p and p r . 
i 

It will be noted that the decomposing constituents are selected to match 
TGA data and have no explicit connection to the identification of actual 
molecule types, of course, it frequently happens that constituents clearly 
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identifiable on TGA traces do represent particular molecules or groups of 
similar molecules, but this is a secondary matter of no significance to the 
in-depth solution. 

4.4 SPECIFIC HEATS 

The transient in-depth solution naturally requires data on the specific 
heats of the virgin material and the char. These are supplied to the program 
as tabular functions of temperature. 

For partially degraded material with P c < p < the analysis of Sec- 

tion 3. 2. 3.4 above has introduced the model that the specific heat is a par- 
ticular hind of mass mean specific heat. In this model, partially degraded 
material is thought of as a mixture of pure char and pure virgin material. 

The specific heat (and hence the enthalpy) of partially pyrolyzed material 
is thus a mass weighted average of the virgin and char specific heats. 

This model is convenient and plausible. Furthermore, the specific heats 
of virgin and char material usually differ only slightly, so that the choice 
of a mixture model is not very important. 

4.5 HEAT OF FORMATION 

Enthalpies for char and virgin material are computed as integrals of the 
specific heat function from the datum temperature to the temperature of inter- 
est, plus a heat of formation. Equations (22) and (23) above indicate this 
calculation. Hence, the user of the computer program must supply heat of for- 
mation data. For most materials such data may be found in various tabulations, 
or may be derived from heat of combustion or heat of pyrolysis data. in the 
case of materials for which no measured data are available, the required heats 
of formation may usually be obtained by theoretical methods or from rules of 
thumb regarding the amount of energy associated with each type of chemical 
bond. For most materials the heats of formation of char and virgin plastic 
have only a minor effect on any ablation results, and hence, great accuracy is 
not required. 

The pyrolysis gas heat of formation is much more significant and will be 
discussed in Section 4.7. 

4.6 HEAT OF PYROLYSIS 

Inspection of the energy Equation (39) reveals that the "heat of pyrolysis" 

or energy effect associated with the degradation of virgin material to char 

is given by (h - h ) ,Btu per pound of pyrolysis gas formed. The computer 
*3 


36 



program automatically calculates this temperature dependent quantity as needed, 
employing the input specific heat functions and heats of formation. 

4.7 PYROLYSIS GAS ENTHALPY 

The in-depth energy Equations (38) and (39) indicate that it is most con- 
venient to have the pyrolysis gas thermal properties as a simple tabulation 
of enthalpy versus temperature. If it can be assumed that the pyrolysis gases 
are in chemical equilibrium as they pass through the char, it is a simple matter 
to run a series of equilibrium chemistry solutions with any convenient compu- 
tation scheme to obtain the necessary enthalpy table, provided that the ele- 
mental make-up of the gas is known. The required elemental data can be ob- 
tained from virgin and char density and elemental analysis measurements. 

Unfortunately, such equilibrium calculations often indicate that solid 
carbon particles should precipitate out of the pyrolysis gas. Allowing these 
particles to collect on the char structure complicates the analysis enormously 
(see the sixth report of the present series) . Thus the analyst is tempted to 
either 

(1) Allow the solid carbon to form, but require it to travel out of 
the char with the gases 

(2) Suppress carbon formation in the enthalpy calculation, thus intro- 
ducing a kinetic consideration. 

Both of these techniques are often tried in practice. 

To do an accurate calculation of the pyrolysis gas enthalpy accounting 
for all possible kinetically controlled reactions would be very difficult 
and would require information about the initial molecular configurations in 
the pyrolysis gas, information not usually available. It is important to 
recognize, however, that the in-depth solution only requires an enthalpy- 
temperature function, so that any physical model for which data are available 
can be accommodated without any changes. in this sense the treatment of the 
pyrolysis gas kinetics is quite general, excluding only carbon deposition 
(" coking") and char erosion. Changes of gas enthalpy with pressure are sec- 
ondary and are ignored. 

It is unfortunate that energy absorption by the pyrolysis gases as they 
pass through the char is very important for most materials of interest. 

Hence, it is important to choose a physical model which is accurate in obtain- 
ing the gas enthalpy curve . 
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4.8 THERMAL CONDUCTIVITY 


Thermal conductivity, as a function of temperature, is required for the 
virgin material and the char. Conceptually, this presents no difficulties 
although in practice good data are very hard to obtain. (Many of the applica- 
tions of the computer program are for backing-out conductivity data from 
measured thermocouple responses taken from instrumented test samples.) 

The treatment of conductivity for partially degraded material is a con- 
troversial, interesting, obscure, and very important topic. In the past it 
it was deemed sufficient to follow the model used for specific heat and to 
use a linear interpolation on virgin material mass fraction; 

k = xk + (1 - x)k 
P c 

where the partially degraded material is regarded as a mixture of virgin mate- 
rial and pure char. The quantity x is the pounds of virgin per pound of 
partially degraded substance. It equals unity for virgin plastic and zero 
for pure char, and at intermediate densities can be shown to be 


x 



1 


°c 

p 


Recent experiments have shown, however, that while this model may be a 
good one for specific heat it can sometimes be a very poor one for thermal 
conductivity. Many materials display, shortly after decomposition begins, 
a conductivity appreciably lower than either the virgin conductivity or the 
char conductivity, a circumstance which cannot be predicted by the model 
described above. Therefore, a recent modification to the computer program* 
has introduced a more general expression 


k = f. (x) k + f (x) k 

where f and f g can be input as general tabular functions of x. With this 
device, formerly inexplicable in-depth thermocouple response could be "pre- 
dicted” quite accurately. It appears that the surface temperature history is 


*This work is reported in Rindal, R. A., Clark, K. J., and Moyer, C. B.: 
Experimental and Theoretical Analysis of Ablative Material Response in a 
Liquid Propellant Rocket Engine, Fourth Quarterly Progress Report, NASA 
Lewis Research Center, Cleveland, Ohio, Contract No. NAS3-7945, 1966. 
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little affected by the new model for mixture conductivity, but the temperature 
profile through the narrow pyrolysis zone is drastically altered for those 
materials in which the conductivity depression is marked, in particular the 
rapid decomposers such as nylon phenolic. 

4.9 SURFACE EMISSIVITY e 

The role and importance of this temperature dependent quantity need no 
elaboration. The chief difficulty is to obtain good data. If partially de- 
graded material is at the surface, the emissivity of virgin material and char 
may be averaged in the same manner as the specific heat. 

4.10 HEAT OF ABLATION, HEAT OF COMBUSTION 

The reader will note that no such quantity as heat of ablation or heat 
of combustion is required as input information for the surface energy balance 
options of the analysis. It will be recalled that in these options the in- 
depth solution is coupled to a general chemistry solution. Heats of ablation 
and heats of combustion are empirical devices designed to circumvent the chem- 
istry solution. These approaches have their advantages, but they are diffi- 
cult to generalize. The only reliable general procedure involves a complete 
surface chemistry solution. With such a solution in hand it is possible to 
calculate the heat of ablation or heat of combustion as a matter of interest 
and for correlating purposes. 

4.11 SURFACE SPECIES 

The identity of the chemical species existing at the ablating surface 
is usually not evident a priori for complex materials. This information is 
not required as input, however. It is a consequence of the chemistry solution. 

4.12 CONCLUSION 

The analysis has been formulated to require the input of fundamental 
physical information in such a way as to preclude any inconsistencies in 
energy terms (such as enthalpies) and in mass terms. Within these restric- 
tions, only input of directly measureable data is required. The three items 
which are probably most important in their influence on computed results, 
the gas enthalpy, the thermal conductivity, and the surface emissivity, are 
unfortunately the items most difficult to obtain good data for, but this, 
of course, is an inherent feature of the problem, not a peculiarity of the 
analysis presented. 
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SECTION 5 


NOTES ON FILM COEFFICIENT MODEL OF THE BOUNDARY LAYER 
WITH HEAT TRANSFER, MASS TRANSFER, AND CHEMICAL REACTION 

5.1 INTRODUCTION 

Section 3.4 above described in general terms how the in-depth energy solu- 
tion can be coupled to the heated surface energy balance boundary condition. 

This coupling is made, through the energy balance, to some computational scheme 
which accounts for the basic transport physics (momentum, energy, chemical 
species) and the basic conservation laws for the boundary layer. 

The boundary layer integral matrix procedure developed under the present 
contract and presented in coded form as the BLIMP program consitutes a very 
accurate and complete procedure for generating boundary layer transport solu- 
tions suitable for this coupling. As might be expected, however, this complex 
routine consumes much valuable computer time for each solution. Hence, economic 
considerations encourage the use of some simpler boundary layer procedure when- 
every possible. 

For this purpose the basic in-depth charring material solution procedure 
has been coupled to a film coefficient boundary layer model which offers speedy 
and economical approximate solutions which still retain all of the essential 
chemical features of ablation events. Thus the in-depth program exists in 
coupled form in two distinct versions, one coupled to the BLIMP program, and 
one coupled to a film coefficient model. Section 3.4 above describes the cou- 
pling in general terms applicable to both coupled versions. For details on 
the BLIMP program and the associated integral analysis the reader may refer 
to the summary report of the present series (in preparation - see Foreword) 
and to Reference 35. Unfortunately, a corresponding detailed presentation of 
the fundamentals of the film coefficient model used for the other coupled ver- 
sion has not previously been presented. Therefore the present section offers 
a few remarks on the film coefficient model as supplementary information. 

Since the chief purpose of the present report is to describe the in-depth solu- 
tion procedure, the following remarks on the film coefficient surface model 
will be in outline form, with arguments and developments omitted. 

5.2 GENERAL POSITION, IMPORTANCE, AND HISTORY OF FILM COEFFICIENT MODELS 

APPLIED TO MASS TRANSFER WITH CHEMICAL REACTIONS 

5.2.1 General Remarks 

Probably the great majority of boundary layer heat and mass transfer 
solutions with chemical reactions have been built on either 
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(1) An empirical approach for any given problem, or 

(2) A complete, accurate solution applied to specific, restricted 
chemical systems. 

Each of these alternatives can be well defended for particular problems. In 
the development of large and costly computer programs, however, it is important 
to build in as much generality as possible. It has been known for a long time 
that film coefficient boundary layer models, originally developed for heat 
transfer calculations, can provide much generality by an extension to include 
mass transfer effects. The film coefficient approach provides the necessary 
transport information, and some appropriate chemical solution routine can pro- 
vide the necessary chemistry solution. 

The fundamental analysis involved in using such a film coefficient scheme 
involves two key parts: 

(1) Derive the film coefficient mass transport relations, and 

(2) Derive the film coefficient energy transport relation 

The actual use of the formulation requires another step, that of calculating 
the film coefficients Pe U e^H an< ^ ^e U e^M" This final aspect is not simple, 
especially for high rates of mass transfer (blowing) from the surface, and 
indeed for strongly non-similar boundary layers probably cannot be carried out 
in a useful way. Nevertheless this task can be carried out for a great many 
problems of importance. 

The first two tasks cannot be carried out without specific consideration 
of the calculation of the film coefficients unless rather loose, heuristic 
arguments are accepted. The number of cases in which Steps (1) and (2) can 
be carried out with precision is rather small. 

In one particular case, however. Steps (1) and (2) can be carried out: 
if Le = 1 and all mass diffusion coefficients are equal, then for a great 
many problems it is possible to calculate p^u^C^, demonstrate that C M “ C H' 
and carry out Steps (1) and (2) without difficulty, regardless of the complexity 
of the chemical events involved. References 36 and 37 provide many examples 
of combined heat and mass transfer for this case, covering a very wide range 
of application areas. This simple model is included as part of the coupled 
computational procedures, but since it appears as a special case of a more 
general formulation to be described below, any further discussion will be post- 
poned . 

Although the Le = 1 and = 1 case provides very powerful sim- 

plicities in analysis, it is quite frequently not an accurate representation 
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of the physical situation in rocket nozzle and reentry problems. The very 
wide range of molecular weights of the chemical species involved makes the 
equal diffusion coefficient assumption inaccurate. 

Therefore, there is a strong incentive to extend the film coefficient 

formulation to cases in which, as a first improvement, C__ / C T1 , and as a 

M H 

further improvement, the mass diffusion coefficients are not equal. In the 
sections below we shall consider these improvements in turn. 

5.2.2 History of Film Coefficient Model for C.. / C 
M H 

It can be shown, for the case in which Le ^ 1 but all the mass diffusion 
coefficients are equal, that it is generally possible and useful to write the 
diffusive flux of a chemical element away from the wall to the boundary layer 
as 


\ = ?e u e C M ( \ - \ 1 < 63 > 

w we 

Extended discussions of this relation may be found in References 24, 26, 37, 
and 38 . 

The corresponding energy transfer relation is not so easily established, 
unfortunately. Exact arguments can be carried through only for the case 
Cm ~ Cjj? for the present case some limiting or perturbation type arguments 
can be carried through for special problems. For example, the case of no net 
mass transfer and a frozen boundary layer is considered in a heuristic way 
in References 37 and 39. Spalding gives a summary of the state-of-the-art 
for this special problem in Reference 40, along with an extensive bibliography. 
Lees discusses this case in Reference 41, along with a number of other cases 
to be mentioned later, as do Fay and Riddell for the specific case of air in 
Reference 42. For "almost frozen” boundary layers (still with no net mass 
addition) Spalding gives brief discussions of the recommended film coefficient 
energy equation in References 37 and 40. For very reactive boundary layers 
with no net mass addition, discussions are given by Lees in Reference 41, 
quoting only Fay and Riddell's results for air in Reference 42. 

For the more important case of net mass addition to the boundary layer 
with chemical reactions, the discussions presented in the literature for the 
case C M t C H appear to be limited to Lees in Reference 41 for the frozen 
boundary layer and to References (26) and (38) . The Lees analysis considers 
the governing boundary layer differential equations and derives, for the frozen 
boundary layer case, an expression for the rate of energy transfer to the wall 
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in which heat and mass transfer coefficients have been identified. The develop- 
ment of References (26) and (38), most completely presented in Reference (38), 
features a looser argument based on the similar appearance of the momentum, 
mass, and energy differential equations of the boundary layer. The analysis 
of Lees and the argument of Reference (38) lead to essentially the same sur- 
face energy equation, which will be discussed in Section 5.3 below. 

5.2.3 History of Extension to Unequal Mass Diffusion Coefficients 

Reference 26 proposes a film coefficient model extended to account for 
the effects of unequal diffusion coefficients. This extension, which appears 
to be unique, is discussed in Section 5.3 below. 

5.3 DISCUSSION OF FILM COEFFICIENT EXPRESSIONS 


5.3.1 Mass Transfer 

For equal mass diffusion coefficients, the diffusive mass transport rate 
of an'element k away from the surface is given by Equation (63). This equa- 
tion can be "derived” for many cases of interest by a solution of the govern- 
ing boundary layer differential equation. In such a derivation, C^ takes 
on the meaning of a convenient collection of physical quantities which is 
essentially independent of the "driving forces" ?^k w a nd . Of course. 

Equation (63) can define or correlate C^ in all problems, but it is generally 
preferable to adhere to the first interpretation. 

For unequal mass diffusion coefficients, Reference (26) shows that if the 
mass diffusion coefficients are related to each other in a particular way 
(which, in fact, appears to represent reality very closely) then it can be 
hypothesized from a study of the relevant differential equations that the dif- 
fusion rate of a chemical element at the wall should be represented by 


where 


X = <>e u e C M<^ - X > 


w e 


\ 4 


Z 


(64) 


(65) 
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(6b) 


and 
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since ft^x^ — 


(67) 



( 68 ) 


The factors F ^ in these equations derive from the particular relation between 
the binary diffusion coefficients which must hold if the governing differen- 
tial equations are to reduce to the forms from which Relation (64) can be 
inferred. This relation is 


ij ' F.F. 

i 3 


(69) 


and can be regarded as an accurate correlation of experimental data for the 
binary diffusion coefficient *^ij* T ^ e < 3 ua ntity D is a constant for a 
given pressure. The constants depend weakly on temperature. 

Thus the rate of diffusion of a chemical element away from the wall is 
given by either Equation (63) if the diffusion coefficients are equal, or by 
Equation (64) for unequal diffusion coefficients. Note that for equal diffu- 
sion coefficients so that Equation (64) reduces to Equation (63) , as 

it must. 

Equations (63) and (64) thus provide the specification of the diffusive 
fluxes necessary to the surface mass balance operations. The overall surface 
mass balance for an element k thus becomes for an ablating, pyrolyzing mate- 
rial (in the absence of condensed phase removal) 



+ 



- % > + ( P v) w\ 

we w 


(70) 


Use of this equation will be discussed in Section 5.4 below. 

5.3.2 Energy Equations 

The expression for the diffusive energy flux from the boundary layer to 
the wall obtained by Reference 26 for equal diffusion coefficient (which is 
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very similar to the less general expressions derived by Lees) is 


q diff 


P U C TT 

e e H 


T b Jb 

H sr - h sw 


Z T, 

(K. - K . ) h . 

1 ' 


T 

w 


(71) 


In this equation it will be noted that the enthalpies of the edge gas and 
the surface gas have been divided into "sensible" and "chemical" parts. This 
division has been necessary in order to assemble the boundary layer differen- 
tial equation terms into groups which resemble groups in the mass and momentum 
differential equations. This resemblance leads to Equation (71) as an energy 
model of Equation (63) . However, in order to establish the resemblance it 
turned out to be necessary to say that the split between "chemical" and 
"sensible" enthalpy occurs at T w - All enthalpy at this temperature is to be 
called "chemical" enthalpy and all additional enthalpy at temperatures above 
this "split temperature" or "base temperature" is the "sensible enthalpy" 
contribution. The base temperature is not to be confused with the "datum 
temperature," at which the total enthalpy of certain basis chemicals (usually 
the elements in certain defined states of aggregation) is zero. 

This enthalpy splitting artifice employed to obtain Equation (71) need 
not be retained any longer however. The term H Sr , measured above the base 
temperature T b = T , may be written 


- E K i. 


,T T 

C„ dT + I W C„ dT + h. 

P i J P i 


+ (rf)KE 


- £ *i 


T T 

f w c dT + h.° 

J P i 


( 12 ) 


where the total (sensible + chemical) enthalpy of the edge gases at the wall 
temperature, h , has been added and subtracted in order to form the difference: 
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edge gas 


(73) 


Recalling that 


sr 


0, we may write 



edge gas 


(74) 


where the total enthalpies in the above relation may be evaluated above any 
datum temperature. Note that the enthalpy of the edge gas at wall tempera- 
ture is for the chemical composition frozen at the boundary layer edge condi- 
tions, not the equilibrium enthalpy of edge gas at the wall temperature. 

Hence the diffusive flux cj^ff becomes 


q diff 


p u C TT ( H - h ) , + p u c 

K e e H v r w edge H e e M 

gas 




T 

K. )h. W 


(75) 


The T w in the final term is now the last remnant of the splitting operation, 
T w 

but h^ merely means the enthalpy of molecular species i at the wall tem- 
perature. The quantity (H r - h^) edge gas does n °t depend on the choice of 
the base state since each enthalpy is for the same physical material, namely 
frozen edge gas. 

Now the overall surface energy balance may be written as 


i u C (H 
e e H r 


h ) . + 

w edge 

gas 


P U 

^ e e M 




(K. - K. )h. w 

l l i 

e w 


+ B 1 h + B'h - B'h 
c c g g w 


+ a w q rad Fc7€T w ~ ^cond 0 (76) 

Only one final detail deserves attention. Equations (75) and (76) have 
been derived in Reference 38 with rather loose arguments. Any checks on the 
correctness and applicability of these equations are of importance. One key 
check requires the results computed from Equation (76) to be independent of 
the enthalpy datum state, A lengthy argument can be constructed to show that 
this is the case. (Equation (75) alone does depend on the datum state, of 
course, through the ku W term.) 

As a final point it may be observed that the derivation of Equations (75) 
and (76) given in Reference 26 does not include any restrictions about the 
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location of the chemical reactions. It cannot be firmly concluded, however, 
that Equation (75) holds for a strongly reacting boundary layer (except, of 
course, in cases where Equation (75) is used to define or C^) . Indeed, 
the separation between the frozen term for heat transfer and the enthalpy 
term' for mass transfer suggests that Equation (75) applies only to a frozen 
boundary layer, reactive wall case. Furthermore, frozen boundary layer analy- 
ses lead to equations equivalent to Equation (75) in References 40 and 41. 
Lastly, physical arguments can be constructed to suggest that for highly re- 
active boundary layers it should be impossible to write general equations 
of the form of Equation (75) . The discussion of pages 263-275 of Reference 
37 is of interest in this regard. 


It may be expected that Equations (75) and (76) may apply with good 
accuracy to frozen and nearly frozen boundary layer problems, but are probably 
less accurate for reactive boundary layers. An indication of the accuracy 
of equations such as (75) for highly reactive boundary layers may be obtained 
by referring to the numerical boundary layer solutions of Fay and Riddell 
(Reference 42). Fay and Riddell's results differed only slightly whether the 
boundary layer was treated as frozen with a catalytic wall or if equilibrium 
was assumed throughout the boundary layer. Even though their solutions are 
limited to air stagnation flow it is believed reasonable to generalize their 
conclusion* to include gases other than air and for body locations other than 
the stagnation point. If this generalization is accepted. Equation (75) may 
be employed for equilibrium, as well as frozen boundary layers, with little 
loss in accuracy. 


Reference 26 extends Equations (75) and (76) to the case of unequal mass 
diffusion coefficients by analogy arguments. The resulting equations resemble 
Equation (64) . The diffusional energy flux to the wall becomes 


Miff 


p e u e C H (H r " 
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w' edge 
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(77) 


and the overall energy balance becomes 


*" . . . the heat transfer is almost unaffected by a nonequilibrium state of 
the boundary layer provided the wall is catalytic and the Lewis Number near 
unity. " 
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( 78 ) 


p u C tt (H - h ) , 4* p u C w 
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+ a< 3«j - FoeT 4 -q , = 0 

w^rad w ^-cond 


Equations (77) and (78) are somewhat less well founded than the corres- 
ponding Equations (75) and (76) , but the equations do have plausibility. 
Furthermore, the equations reduce to Equations (75) and (76) for equal dif- 
fusion coefficients, as they should, and Equation (78) can be shown to be 
independent of the enthalpy datum state and thus fulfills a basic physical 
requirement . 


5.4 USE OF FILM COEFFICIENT EQUATIONS TO CALCULATE SURFACE 
ENERGY BALANCE 


5.4.1 Basic Aspects 

The film coefficient model has provided simple expressions for the dif- 
fusive transport rates of mass and energy through the boundary layer to the 
wall. It is now necessary to indicate how these relations may be expeditiously 
used in computing the surface energy balance. Section 3.4.2 above describes 
how it is useful before doing the surface energy balance, to have the diffu- 
sive flux terms and h , q*, and T as functions of m and m before 

w w eg 

forming the energy balance. For this purpose, the film coefficient mass 

transport relations may be coupled to some convenient chemistry routine which 

solves the coupled mass balance + chemistry problem. Equations (75) and (78) 

su 99 es t that m and rn can be normalized on p u C w to increase the qen- 
c g s e e M 3 

erality of this procedure. and B^ can then be regarded as independent 

variables in the construction of a table of values of T , Y zt h. w , and 
Z. h. for given pressure, edge composition, char composition, and gas 
composition. (The term q* may be lumped with this last term if desired.) 

The resulting table meets all the requirements outlined in Section 3.4.2 

and can be used for an iterative determination of that B’ which satisfies 

c 

the surface energy balance, given from the in-depth solution. 

The Aerotherm Equilibrium Surface Thermochemistry Program (EST) is a 
coupled mass balance - equilibrium chemistry routine which generates such 
tables for chemical equilibrium conditions. Reference 43 describes this 
program. A somewhat newer version with improved computational techniques is 
denoted the Aerotherm Chemical Equilibrium Program (ACE) and is described 
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in Reference 44. This program allows certain reactions to be kinetically con- 
trolled and hence has more generality. 

5.4.2 Input and Correction of Heat Transfer Coefficient 

To employ the film coefficient formulation just described, the program 
user must provide the program with values of the heat transfer coefficient 
p^u^C^ as functions of time. Two practical problems must be settled in this 
respect: 

(1) How is c related to C„? 

M H 

(2) Can both and be specified as functions of edge conditions 

(i.e., of time) independent of the subsequent problem solution 
(i.e., mass transfer rates and body shape)? 

In answer to the first question it may be stated that within the present 
formulation it is adequate to take the ratio C„/C_, as a constant. The value 
of this constant is a measure of the ratio of the mean mass transfer aspects 
of the boundary layer to the mean heat transfer aspects. For equal mass dif- 
fusion coef ficients, a vast amount of experimental data (as summarized, for 
example, in Reference (37)) suggest the correlation C-./C^ = Le Y . It may be 
hypothesized that for unequal mass diffusion coefficients the same procedure 
may be employed with the Lewis Number, Le, defined by the procedure set forth 
in Reference 26 involving D. Thus, the input to the program consists of a 
time table of values for P e u g C M and the constant factor C M /C H . 

The answer to question (2) , changes of C„ with body shape are occasionally 

H 

of interest and may easily be accounted for. A more important problem concerns 
the dependence of C 0 on the actual rate of mass transfer. This problem has 

n 

been ignored up to now, the implication being that is determined by the 

boundary layer edge aerodynamics alone. This is known to be incorrect. The 
value of depends fairly strongly on m. If we denote the C with m = 0 

as c Ho' this dependence is shown by both data and analysis to be accurately 
represented by 


where 
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(79) 
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Equation (79) accurately correlates a large amount of laminar data if X is 
chosen as 0.5. For turbulent flow a X of 0.4 appears to be slightly better. 
The computer program provides for X as an input parameter and then automati- 
cally computes the "blowing reduction." 

5.5 CONCLUSIONS 

A film coefficient model for boundary layer heat and mass transfer with 
chemical reactions can be developed for the governing differential equations 
of the boundary layer. Of the resulting film coefficient expressions, those 
for mass transport can be coupled to any convenient chemistry routine to pro- 
vide boundary conditions for an in-depth response calculation. The boundary 
information and the in-depth calculation are coupled through the film coef- 
ficient surface energy balance, providing a complete link between the boundary 
layer and the sub-surface material. 

The film coefficient model of the boundary layer provides a powerful, 
economical, and quite general alternative to empirical procedures and to 
elaborate "exact" procedures. The model can be proved to be accurate for a 
wide range of problems, and appears to apply with good accuracy for an even 
wider range of problems for which good accuracy proofs cannot be constructed. 


SECTION 6 
PRESSURE DROP 


The pressure drop across the char layer resulting from pyrolysis gas 
flow is of interest because it may give rise to excessive char stress levels. 
Knowledge of the pressure distribution through the char is necessary for evalu- 
ating char layer stress levels. An empirical relation for evaluating the 
pressure distribution through the char layer is rationalized and presented 
first, in Section 6.1 and is followed, in Section 6 . 2 by the finite difference 
formulation . 

6.1 PRESSURE DROP CORRELATION EQUATION 

In order to obtain an expression relating local pressure in the char 
layer to other pertinent variables it is useful to examine experimental data 
tahen for a range of variation of the pertinent parameters of interest, and 
to deduce the most significant effects by generalizing this data. The data 
presented by Green in Reference 45 is employed for this purpose. Green 
presents a compilation of data for flow of various gases through a wide variety 
of porous media for a large range of flow conditions. In order to relate 
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pressure gradients to viscous and inertial forces he employs the correlation 
equation first proposed by Reynolds: 


dP 

dy 


aptv + gpv 3 


(81) 


where the first term represents viscous forces and the second, momentum forces. 
The quantities a and g are empirical coefficients, fj is the gas viscosity, 
and p is the gas density. The velocity (v) in Equation (81) is a "super- 
ficial velocity” defined on the basis of the gas flow rate per unit projected 
area in a plane normal to the velocity vector. 


v 


A 


m 

pA 


(82) 


Referring to the above equations, the ratio of inertial to viscous forces may 
be written as: 


Inertial Force = p g 
Viscous Force ” cx|iA 

The compilation of data presented by Green includes tabulations of the empir- 
ical coefficients a and g for a wide variety of porous media including 
packed beds of irregular and spherical particles ranging in nominal size from 
0.08 inch to 0.1875 inch and for close packed wire screens ranging from 60 
to 5 mesh. The nominal range of porosities included varies from 0.3 to 0.8. 
The experimentally derived coefficients a and g are shown in Figure 3 
along with a straight line fit to the data. 

a = 0.794 x 10 4 g 

where a has units of ft -2 and g has units of ft -1 . The correlation is 
not excellent, but appears appropriate for order-of-magnitude considerations. 
Substituting this relation into the equation obtains 


inertial Force , 26 x 10’ 4 ^2 
Viscous Force ’ (iA 

For the temperature range of interest (2000°R - 5000° R) the gas viscosity will 
range from about 0.3 x 10~ 4 to 0.5 x 10 -4 lb/ft-sec. Considering a gas vis- 
cosity of 0.4 x 10 -4 lb/ft-sec, the ratio of inertial to viscous forces may 
be written: 
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Inertial Force 
Viscous Force 


= 3.15 


m 

A 


where m ^/A has the units of lb/ft 8 -sec. It may therefore be concluded that 
for pyrolysis gas flow rates of 0.01 lb/ft 2 -sec and less, the inertial terms 
may be ignored in which case Equation (81) reduces to Darcy's Law. 


dP 

dy 


djJV 


where a” 1 is the permeability. It is believed appropriate, however, that 
the more complete correlation Equation (81) be employed for pressure drop 
calculations in the char layer of ablating materials since it is valid for 
a wider range of conditions of practical interest. 


6.2 FINITE DIFFERENCE FORMULATION 

It is desired to cast the pressure drop correlation Equation (81) into a 
form including only variables which are readily available from the computation 
scheme for solving the mass and energy conservation equations (Section 3.2.4). 

The velocity of the gas passing through the porous char is related to 

the gas flow rate per unit area by Equation (82) . Introducing the gas state 

equation (P = p RT/r) into Equation (82) yields: 

y y 
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(83) 


It is convenient to define a viscosity law for the pyrolysis gas 


M- 


(84) 


Substituting Equations (83) and (84) in (81) yields the following rela- 
tion for a linear pressure drop across a node (n, see Figure 2) . 


dP 

dy 
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where the coefficients a, 3, and the gas molecular weight, ft? g are input 
as functions of temperature and m g represents the pyrolysis gas flow rate 


entering the node n. Equation (85) is solved for n = 1,2,3 


down to 


the first node in the virgin material where 


%i = °- 


The pressure at the 
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top of node 1 is taken equal to the prescribed boundary layer edge pressure 
(Pi = p^) . The pressure at the top of any other node is related to the pres 
sure drop across and pressure at the top of the node above it. 
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( 86 ) 



i=i 


(87) 


The computation is performed in an explicit manner, that is, after obtaining 
solution of the energy and mass conservation equations, a single pass through 
the nodal network from heated surface (n = 1) to the virgin material yields 
the pressure distribution entirely in terms of conditions existing at the end 
of the' previous time step. 


SECTION 7 

SOME OPERATIONAL DETAILS OF THE AEROTHERM 
CHARRING MATERIAL ABLATION PROGRAM, VERSION 2 

7.1 INTRODUCTION 

While it is not the purpose of the present document to offer a detailed 
description of the charring material computer program which is based upon 
the analysis presented in Section 3 above, it does seem useful to include 
a few words of program description at this time. Complete descriptions, 
user’s manual, and flow charts are given in References 46, 47, and 48. 

7.2 PROGRAM DESCRIPTION 
7.2.1 General Remarks 

The Charring Material Ablation (CMA) program is a coded procedure for 
calculating the in-depth thermal response of a charring, ablating material. 
The solution is obtained through difference equations discussed in Section 3 
above. This section presents a brief over-view of the program operation 
and output. 
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7.2.2 Pro gram Objectives 


The program produces in-depth temperature and density histories, plus 
surface recession rate as a function of time. In addition to this basic out- 
put, the program outputs a number of integrated energy terms and various mate- 
rial property data of interest. Section 7.2.5 below gives a more detailed 
description of the program output. 

7.2.3 Program Capabilities 

The Charring Material Ablation Program is an implicit, finite-difference 
computational procedure for computing the one-dimensional -transient transport 
of thermal energy in a three-dimensional isotropic material which can ablate 
from a front surface and which can decompose in depth. Decomposition reac- 
tions are based on a three-component model. The program permits up to eight 
different backup materials of arbitrary thickness. The back wall of the com- 
posite material may transfer energy by convection and radiation. 

In one program configuration, the ablating surface boundary condition 
may take one of three forms: 

Option 1 - Film coefficient model convection- radiation heating with 
coupled mass transfer, including the effects of unequal 
heat and mass transfer coefficients (non-unity Lewis 
number) and unequal mass diffusion coefficients. Surface 
thermochemistry computations presume chemical equilibrium 
at the surface. 

Option 2 - Specified surface temperature and surface recession rate 

Option 3 - Specified radiation view factor and incident radiation 
flux, as functions of time, for a stationary surface. 

Any combination of the first three options may be used for a single 
computation. Option 3 is appropriate to cooldown after termination of con- 
vective heat input and is often useful in conjunction with Option 1 and 2. 

In another configuration, the program may be coupled to the Boundary 
Layer Integral Matrix Procedure. In this arrangement, the total assembly 
is designated the CABLE program and is described in other reports of this 
series (see Foreword) . 

The program permits the specification of a number of geometries: 

(1) Plane 

(2) Cylindrical or annular, with heated surface either inner 
or outer 
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Spherical or spherical shell, with heated surface either inner 
or outer 

(4) General "thermal stream tube" geometry, area varying as depth 
to any power 

(5) General "thermal stream tube" geometry, area varying arbitrarily 
with depth. 

The rear surface of the last node may be specified as insulated, or may 
experience convective and radiative heat transfer to a "reservoir" at a spe- 
cified reservoir temperature if a rear surface convection coefficient and an 
emissivity are input. 

Material properties such as thermal conductivity, specific heat, and 
emissivity are input as functions of temperature for virgin plastic and char. 
For partially decomposed material, the program perforins an appropriate averag- 
ing to determine effective material properties. 

7.2.4 Solution Procedure 

The basic solution procedure is by a finite difference approach. For 
each time step, the decomposition relations are solved and then the in-depth 
energy fluxes constructed in general terms. These are then harmonized with a 
surface energy balance (if a surface energy balance option is being used) and 
the in-depth temperatures determined. New material property values are set 
up and the solution is ready for the next time increment. 

7.2.5 Output Information 

The CMA program outputs instantaneous mass ablation rates and blowing 
parameters for char and pyrolysis gas, total integrated mass ablation of char 
and pyrolysis gas, total recession and recession rates of surface, of the char 
line, and of the pyrolysis line. It also outputs the surface energy flux 
terms, namely, the energy convected in, energy radiated in, energy reradiated 
out, chemical generation, and conduction away (q conc j) • Further, it describes 
how the input energy of < l conc 3 i s "accommodated" or "partitioned" in the 
solid material. Part of the energy is consumed in decomposing the plastic, 
part is consumed in sensible enthalpy changes of the solid, and part is 
"picked-up" by the pyrolysis gases as they pass through the char. 

Finally, nodal temperature, density and enthalpy are given for every de- 
sired output time, as well as a description of the current material state for 
each node. 
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Thermocouple and isotherm output can also be called for. Figure 4 shows 
a representative page of output from a problem solution. 

7.2.6 Operational Details 

7. 2. 6.1 Storage Requirements 

The storage requirements for the CMA program depend strongly upon the 
coupling mode in use. Coupling to a film coefficient model for the surface 
energy balance involves so much table storage that the program will barely 
fit a 32 , 000 word machine with full table sizes. In certain cases a reduc- 
tion in table sizes will allow the program to fit on a smaller machine. 

Coupled to the CABLE program, which eliminates the need for storing extensive 
boundary condition tables, the CMA program requires less than 8000 words of 
storage . 

7 . 2 . 6 . 2 Running Time 

Computation time depends, of course, on the problem being computed, but 
experience to date indicates that CMA computations run in roughly 1/3 of real 
time for 11 typical" charring material problems, for machines of the IBM 7094 
speed class. 

7.3 SAMPLE PROBLEM SOLUTIONS 

7.3.1 Some Typical Problems 

As an illustration of the general performance of the charring material 
computer program, Figure 5 presents a graphic representation of the in-depth 
density history for a nylon-phenolic heat shield material exposed to a typi- 
cal air reentry environment. The peak cold wave heat flux was about 800 
Btu/ft s sec. Figure 6 shows some in-depth thermocouple temperature response 
predictions compared to test data. Figure 7, from a different problem, 
shows a machine made plot generated by a plot routine coupled to the CMA 
program. 

7.3.2 Additional Examples 

Additional examples of the program performance may be found in Appendices 
B and C. Appendix B, which is an early technical note written during the de- 
velopment of the program, reports the results of heat conduction solution 
tests. Appendix C describes some computational experiments exploring interest- 
ing aspects of solution smoothness for various charring material problems. 
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SECTION 8 

SUMMARY AND CONCLUSION 


8 . 1 GENERAL REMARKS 

The preceding sections have described the analysis and an associated 
computer program for the calculation of the in-depth response of a charring 
jl pyrolyzing material. The general objective of the development effort has 
been to produce a computation scheme which accounts for those physical events 
common to a wide range of technically important applications, so that the re- 
sulting program has as much generality and flexibility as possible. To this 
end, the analysis accounts for the basic in-depth pyrolysis problem, which 
is observed to be common to a wide range of problems, and excludes coking 
(char densif ication) , subsurface char erosion by pyrolysis gases, thermal 
expansion, condensed phase char reactions, and mechanical damage mechanisms. 
All of these are particular to particular materials, or material types. For 
such materials, the basic program can be modified to include these special 
e f fects . 

The basic program generates a one-dimensional in-depth solution, but the 
cross sectional area of the material analyzed may vary with depth (thermal 
stream tube) . pyrolysis may occur through three distinct Arrhenius-type 
kinetic reactions. 

An important feature of the program is the range of physically realistic 
boundary conditions available for the heated surface. These include 

(1) Specified temperature and recession rate 

(2) Radiation energy balance with zero recession and no convection 
(cool down or soak out) 

(3) Coupling through a film coefficient model to surface thermo- 
chemistry solution including general heterogeneous equilibrium, 
or heterogeneous equilibrium modified by certain rate controlled 
reactions, both models including the effects of the melting and 
total removal of surface species formed at temperature above 
their melt or fail temperatures 

(4) Coupling to a general, non-similar, reacting boundary layer 
solution including homogeneous and hetergeneous kinetic effects, 
with surface melting or failing. 

This range of possibilities offers opportunities for economy during routine 
in-depth studies or during computations for which film coefficient models 
are adequate, while preserving the capability of doing very accurate coupled, 
simultaneous boundary layer and in-depth solutions. 
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8.2 EXPERIENCE WITH THE IN-DEPTH SOLUTION ROUTINE (CMA) 

The in-depth solution routine has been applied to a fairly wide variety 
of materials and some brief account of these applications may provide some 
useful orientation. By suppressing pyrolysis effects, the program can be 
used for simple transient heat conduction problems? applications here have 
included such refractories as alumina, boron nitride, tungsten, and graphite. 
The program has been used very extensively for graphite and carbon phenolic 
rocket nozzle studies. These computations include numerous series of speci- 
fied temperature and recession rate runs to back out thermal conductivity 
data for these materials, as well as a very large number of coupled solutions. 

The program has been run on many parametric studies of nylon phenolic 
with generally excellent results. The very rapid decomposition rate of this 
material occasionally has caused oscillations in the pyrolysis gas flow rate, 
with attendant minor oscillations in surface temperature and surface reces- 
sion rate. Fortunately, overall aspects of the solution, such as total reces- 
sion, mean temperatures, and isotherm penetration seem unaffected by these 
oscillations. Appendix C gives an extensive discussion of the oscillation 
problem. In one pathological case for nylon phenolic, a step application of 
a very high heat flux caused a total disintegration of the solution process. 
This disintegration appeared to be due to the explicit decomposition treat- 
ment combined with the implicit energy treatment and could be circumvented 
by a slight softening of the step transient. No other solution breakdowns 
have been observed with the program. 

The program has also been applied to asbestos phenolic with good results. 

The program has been extensively used for silica-containing insulators 
such as silica phenolic. The results here have been excellent so long as the 
main assumptions of the analysis are not violated. Materials with much silica 
occasionally display physical events such as thick liquid layer runoff and sub- 
surface char reactions, for example condensed phase silica-carbon reactions. 
Since the program was not designed to account for these events it is not al- 
ways applicable to these materials. 

The program was once used to predict the freezing of water, for which 
case the freezing reaction was input as one of the Arrhenius reactions avail- 
able in the program structure. The results were quite accurate . 


8.3 CONCLUSION 

The analysis presented here and programmed as the CMA program has been 
applied to a wide range of materials of technical interest. The results 
have been excellent and the program appears to be thoroughly checked out and 
operational . 
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Figure 2. Finite difference representation. 
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0 Total Surface Recession * 0.41 inch 
□ Total Char Penetration * 0.68 inch 
A Thermally Unaffected Zone Below 0.79 inch 
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Figure 6. Predicted Temperature Histories at Surface and 
Selected Thermocouple Locations, Nylon Phenolic 
Re-entry Problem. 
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Figure 7. Typical Machine Plotted Isotherm Depths 
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EQUATIONS FOR COEFFICIENTS A fi , B n# C , and D 
IN IN-DEPTH ENERGY EQUATION ARRAY 

The coefficients A^, , C n , and in the array of Equation (58) are 

determined by Equations (54), (56), and (57). 

For nodes in the ablating material except the first and last, from 
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For the surface node of the ablating material, to which Equation (54) applies 

A 1 = °' B 1 and c l re tain the forms of Equations (A-2) and (A-3) above, 

while D, has the form of Equation (A-4) plus the term q A /6 
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For the last node of ablating material, to which Equation (57) applies. 
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Here NBM refers to the first node of the hack-up material. if there is 
no back-up material, the resistance 6 NBM^NBM^NBM re P laGec ^ by the back 

wall resistance V h total A^ L where ^ total is an effective heat transfer 
coefficient for convection and radiation. Since T' . will then refer to a 

NBM 

known "reservoir" temperature, the term c nl t nbm may be added to d nl , leaving 
only two unknown temperatures. 

For back-up materials the coefficients are simply 
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6 /2 

n' 

k A 
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T (pC ) 
n 'r p 7 n 
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For the first back-up node, n = NBM and n-1 - NL. For the last back- 
up node (n = N) , n+1 refers to the back-wall condition, described above. 

The term 6 n+1 /k n+1 A n+1 becomes ^/^total^Sl and D N ;i ' S moc *i f: *- ed as described 

above for the case where the last node of ablating material is the last node. 
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GENERAL 

A 

A 

B 


c 

P 

d 

E„ , E 
l * p 

erf ( ) 
er fc { ) 
h 


i 

k 


n 



x 

X* 


area 

dimensionless constant employed in equations (4) and (5) 
constant defined in equation (A-18) 
specific heat 
nodal thickness 

defined by equations (A-7) , (A-8) respectively 

error function 

Q 1 - erf ( ) 

enthalpy 

index 

thermal conductivity 

dimensionless constant employed in equations (4) and (5) 

heat flux per unit area at surface 

radius coordinate (origin at center) 

external radius 

surface recession rate 

solution to difference equation for temperature 
solution to differential equation for temperature 

- V 

distance coordinate with origin tied to receding surface 
similarity variable x/2 '\f~aQ 


SUBSCRIPTS 


c value at center 

e value at edge or surface 

min minimum value in region of interest 

n node index 

0 Initial value, value as x oo 


v 



GREEK 

a 

3 

A<3 

6 

e 

Q 

SPECIAL 

A 


thermal diffusivity k/p 

arithmetic error 

time increment 

slab thickness 

error T* - t* 

time 

equal to by definition 
approximately equal to 



CONDUCTION SOLUTION CHECK-OUT 
OF THE CHARRING MATERIAL ABLATION PROGRAM 


B.l INTRODUCTION 

This report presents the results of a series of check-out computations 
testing the heat conduction aspects of the Aerotherm Charring Material Abla- 
tion (CMA) Program, Version 2. Calculations included transient conduction 
cases for slabs, semi-infinite solids, spheres, cylinders, and receding semi- 
infinite solids. 


B.2 ASPECTS OF PROGRAM 

The program treats the heat conduction and solid convection in a fully 
implicit manner. Thus the difference equation at each finite-difference 
"node" is, for a constant density material. 


pc (T 1 - T ) 
y p n n 


A G 

A d 
n n 


T 1 - T ' 

tw n 

W* ■ V 2 


tjt 1 _ m * 

n n+ i 


k A 
n-i n-i 


k a 
n n 


V 2 

k A 
n n 


d /2 
n+r 

^n+ i A n + i 


+ [ (ph) n +1 + - W " <P h >n " (P c p>n (T A " T n>] 

Here primes denote "new" temperatures at the end of the time step A 9 . 
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B . 3 CONSTANT PROPERTIES SLAB 

A single problem was computed for a 10-inch slab with insulated back 
face, uniform initial temperature, and a step change in front surface tem- 
perature. Property values were constant at 


k = 1 Btu/sec ft °F 
0 ^ = 1 Btu/lb °F 
p = 1 lb/ft 3 . 


The nodal 


distribution 

was , 

from 

the 

surface 

* 

1 

node 

of 

0.10 

inches 

thickness 

2 

nodes 

of 

0.20 

inches 

thickness 

2 

nodes 

of 

0.25 

inches 

thickness 

4 

nodes 

of 

0.50 

inches 

thickness 

5 

nodes 

of 

1.00 

inches 

thickness 

1 

node 

of # 

2.00 

inches 

thickness 
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The maximum time step permitted was 0.002 seconds. 

Figure B-l shows the computed results compared to exact curve presented 
by Schneider (Reference B-l) which have been computed according to the 
solution 



00 



i»i 


(iTT/2) 2 (ote/6 2 ) ilr 

s in "tt x i * 1,3,5.... 


( B-2) 


The agreement between the exact curves and the points produced by the 
program is very satisfactory. 


B . 4 CONSTANT PROPERTIES CYLINDER 

A computation exactly similar to the slab computation described above 
was done for a cylinder with an external radius of 10 inches. Property 
values, nodal sizes, and time step limit were as in Section 3. 

Figure 2 displays the results with computed points shown as open symbols. 
In contrast to the slab case, agreement between computed points and the exact 
solution is not good for the deepest nodes, which respond faster in the cyl- 
inder problem than in the slab problem due to their smaller mass. 

The problem was recomputed with the maximufn time step halved to 0.001 
seconds and with the following nodal size distribution. 

1 node of 0.10 inches thickness, 

2 nodes of 0.20 inches thickness, 

2 nodes of 0.25 inches thickness, 

18 nodes of 0.12 inches thickness. 

The results for this calculation are shown in Figure 2 in the solid 
figures. The exact solution curves are given by Schneider (Reference B-l) 
computed from the relation (Reference B-2) . 


t 


t 


e 



1 





i»i 


T [m. - 

o l 1^1 

W 


( B-3) 


where the M i are the roots of the zero order Bessel function J Q . (For 
checking solutions at r/r^ ratios different from those presented by Ref- 
erence B-l and for greater accuracy in checking, a small program for computing 
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Equation (B-3) was generated. For this purpose, Reference B-3 gives conveni- 
ent Bessel Function expansions as 


(1.) (x) 1 3 


J 0 ( X) 


1 - 2 . 2499997 ( -jj ? + 
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- 0.0039444[y) 10 + 


1 . 265620ejyJ * 
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0 . 000310o( 1P 
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(B-4) 


The quoted error p is < 5 x 10 0 , but Aerotherm experience shows that 
machine computations introduce errors of about 10 _? . 


(2.) x ^ 3 J (x) 
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with 


x 1//p f cos 6 
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0.79788456 - 0.00000077^] - 0 . 0055274o(^j 


- 0.00009512[-] 

( + 0.00137237[-| 

l xi 

l xl 

- 0 . 000728951 — ' 

\ x J 

j + 0.00014476(^] 


P < 1.6 x 10 
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< 7 x 10 , 



( B-9) 


Reference B-3 also gives a useful tabulation of the values for use 

in Equation (B-3) . 

The agreement between the exact solution and the computed points is 
excellent. 


B.5 constant properties sphere 

Two computations were made for a 10-inch sphere, the property values, 
time step limits, and nodal sizing corresponding to both cylinder calculations 
described in Section 4 above. 
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Figure b-3 shows some computed results (as plotting symbols) compared 
to an exact solution. The exact solution is given by Reference B-2 as 
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These solutions are conveniently plotted in Reference B-l. A small program 
for producing results of greater accuracy and at general r/r^ points was 
generated to facilitate comparisons. 

Figure B-3 indicates that for a "proper" choice of nodal sizes, the 
computed results match exact solutions. 


B . 6 SEMI -INFINITE SOLID WITH CONSTANT PROPERTIES 

A semi-infinite slab was simulated with a number of nodes large enough 
to insure that the final node showed no temperature response during the 
computation. Property values were taken as 

k = 1 Btu/ sec ft °R 

c = 1 Btu/ lb °R 
P 

p = 1 lb/ f t 2 

The nodal distribution was, from the surface. 


1 

4 
1 
1 
7 

5 

5 

6 


node of 
nodes of 
node of 
node of 
nodes of 
nodes of 
nodes of 
nodes of 


5 inches 
10 inches 
12 inches 
15 inches 
24 inches 
48 inches 
96 inches 
192 inches 


thickness , 
thickness , 
thickness , 
thickness , 
thickness , 
thickness , 
thickness , 
thickness . 


The maximum time step was limited to 5.0 seconds. 

The exact solution to the semi-infinite solid problem with uniform ini- 
tial temperature t and step surface temperature at time 8=0 is 
a similarity solution 
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1 - erf (x*) 


(B-12) 


where 



(B-13) 


Figure B-4 shows the exact similarity profile compared to two of the 
many computed profiles. The circled points are for the computed response 
at 1 second, after only 8 computational steps. Even at this early time the 
computed profile is very close to the exact profile. The points in tri- 
angles are for time 0 = 64 seconds, after 64 computational steps.* By 
this time the computed points correspond to the exact solution to three 
significant figures. 

As a further check for this problem, exact solutions are available for 
surface heat flux rate. Figure 5 shows the surface heat flux rate and inte- 
grated heat flux plotted versus dimensionless time. (The arbitrary length 
r l ^ as keen introduced in order to non— dimensionalize. ) Agreement between 
the exact solutions of Schneider (Reference B-l) and the computed results 
is excellent. 


B . 7 SEMI-INFINITE SOLID WITH VARIABLE THERMAL CONDUCTIVITY 

Halle (Reference B-4) has presented exact solutions for certain varia- 
ble property problems useful for checking the computer programs. Of these 
solutions, the one for the step wall temperature problem with conductivity 
varying with temperature is 
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A erfc(x*) + (1-A) erfc(nx*) 


(B-14) 


where 




A k o A 

and a ■ and k ■ k(t ) , 

o pc o ' o y 


provided that k varies as 
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Ae x * + (1 - A) — e n x * 

n 

. -x*? -n?x*? 

Ae + (1 - A) n e 


( B-l 5) 


♦Although the maximum time step allowed is five seconds, the program applies 
various limiting criteria during rapid transients so that the average time 
step is less than five seconds. 
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The constants A and n may be chosen to give a useful k(t) function, 
which is obtained by cross-solving or cross-plotting between Equations (B-14) 
and (B-15) . 

Two cases were computed for check-out, one with A = 0.5 and n - 2.0 
(corresponding to Case 4 on Figure 2 of Reference B-4) and the second with 
A = 1.85 and n => 2.0. Nodal sizes, maximum time step, specific heat, and 
density were all identical to those in the constant properties computation 
described in Section 6 above. The conductivity k Q was .taken as 1.0 Btu/ 
ft sec °F . Figure B-6 shows the thermal conductivity variations for the two 
cases, with k/k Q plotted against (t - t Q )/( t e ” * Tabular values used in 

the program are circled. The tabular values used are as follows: 


A * 

Case 1 
0.5, n - 2.0 

A - 

Case 2 
1.85, n =■ 

2.0 

o 

It 

o 

°R, t e - 1000°R 

fc o = 0 


1000°R 

(°R) 

k/k o 

t(°R) 


kA 

0 

1.0 

0 


1.000 

100 

0.91 

88.2 


1.000 

200 

0.76 

165.3 


1.009 

300 

0.66 

287.0 


1.036 

400 

0.60 

366.5 


1.066 

500 

0.56 

457.0 


1.109 

600 

0.53 

555.5 


1.201 

700 

0.52 

656.6 


1.340 

800 

0.51 

753.4 


1.575 

900 

0.50 

838.3 


1.989 

1000 

0.50 

905.3 


2.762 



952.1 


4.306 



968.5 


5.566 



981.2 


7.174 



991. 3 


8.783 



1000.0 


9.500 


During the computation the program obtains thermal conductivity values 
by linear interpolation between these points. 

The first case corresponds to a material whose thermal conductivity de- 
creases rapidly but nearly linearly with temperature, as for example, a pure 
metal, while the second case models a material whose conductivity increases 
very rapidly at high temperatures. This case resembles the charring ablator, 
for which the char layer has a much higher conductivity than the virgin material. 
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Figure B-7 shows the exact and the computed results for the two cases, 
along with the constant properties exact solution to illustrate the degree 
of departure of the two cases from the constant properties case* Agreement 
between the computed and the exact solutions is excellent. 


B.8 SEMI-INFINITE RECEDING SOLID 

Check-out of the convection aspects of the computation requires a problem 
with surface recession. An analytical solution is available for the transient 
response of a semi-infinite solid initially at uniform temperature exposed to 
a step in surface temperature and to a step in surface recession rate S . 

For the constant properties problem it can readily be shown that the tem- 
perature profile approaches a quasi-steady form 


t - 
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t 

0 



(B-16) 


where the x coordinate origin is tied to the receding surface. The work of 
Reference B-5 shows that a useful measure of the approach to steady state is 
provided by the variable 


P c p S(t e ~ t o ) 


comparing the amount of solid convection pick-up to the amount of energy con- 
ducted into the solid. This term is initially zero and approaches unity in 
the steady state. 

Figure B-8 shows the exact transient response, calculated from results 
of Reference B-5, compared to computed results for a problem with the same 
nodal size distribution as the infinite slab problem described in Section 6 
above. The property values were 

k = 1 Btu/ft sec °F 

c = 0.1 Btu/lb °F 

P 

p = 1 lb/ ft 3 

The specified surface recession rate S was 0.5 ft/sec., which limited the 
time step to a maximum value of 1/6 second. (The program will not allow 
more than 20 percent of the thickness of the smallest node to be ablated away 
in a single computational step.) 
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The agreement between computed results and the exact solution is excellent. 

Three separate problems were calculated in order to judge the accuracy of 
the computation of the steady state temperature profile. Parameter assignments 
for these calculations were as follows: 



Case 1 

Case 2 

Case 3 

k (Btu/ft sec °F) 

1.0 

1.0 

1.0 

Cp (Btu/lb °F) 

1.0 

0.1 

0.1 

p (lb/ft 3 ) 

1.0 

1.0 

1.0 

ct (ft 2 /sec) 

1.0 

10.0 

10.0 

S (ft/sec) 

0.5 

0.5 

0.5 

A9 max (sec) 

0.167 

1.6 

1.6 

"Steady State" time (sec) 

100.0 

300.0 

300.0 

Dimensionless "Steady 
State" time S 2 9/a 

5.0 

2.739 

2.739 

Parameter Sd . /a 
mm 

0.5 

0.2 

0.1 


Case 1 

Case 2 

Case 3 

Nodes 

1 at 5" 

9 at 48" 

25 at 24" 


4 at 10" 

5 at 96" 

9 at 48" 


1 at 12" 

6 at 192" 

5 at 96" 


1 at 15" 

5 at 500" 

6 at 192" 


7 at 24" 


5 at 500" 


5 at 48" 




5 at 96" 




6 at 192" 



Figure B-9 shows the results of 

these three 

calculations 

and indicates 


that Case 1 gives a poor steady state profile while Case 2 is improved but 
still not good, and Case 3 gives a good result. The key difference between 
the three cases lies in the relative sizes of a single heat conduction term 
and a single solid convection term in the finite difference energy balance. 
Equation (B-l) . 
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conduction 

term « 

S (T n+ i - V 

(B-17) 

convection 

term sc 

P C P * (T n +1 - V 

(B-18) 

convection 

term 

p y d sa 

( B-19) 

conduct ion 

term 

k a 


The value of this parameter for d . is given in the table above. 

min 

Thus Figure B-9 suggests that an accurate calculation requires 


Sd . 


mm 

a 


< 


0.1 . 


(B-20) 


It cap be shown through an analysis of truncation errors (see Sub-Appendix B--1) 
that this conclusion is not generally valid. it holds only if 

1. The nodes are small enough that truncation errors in the 
conduction terms are "small," 

2. The temperature profile is "close" to e -SX//a . 

Even though the criterion (B-20 1 ) is not generally valid, it is still a 
useful one since many ablation problems approximately satisfy the two conditions, 

B . 9 SUMMARY AND CONCLUSIONS 

The purpose of the test calculations was to show that the Aerotherm 
Charring Material Ablation (CMA) Program, Version 2 (implicit) performs heat 
conduction calculations properly. Checks included computations for slab, 
cylinder, spherical, and semi-infinite solid geometries exposed to step sur- 
face temperature increases and several calculations of semi-infinite solid 
geometries exposed to simultaneous steps in surface temperatures and surface 
recession rate. 

Results indicated that when nodes were sized properly, the results gen- 
erated by the program matched exact solutions and that, therefore, the program 
is judged to operate properly in its heat conduction and solid convection 
aspects . 

Poor choices of node sizes can of course lead to unsatisfactory results, 
but a complete analysis of this subject is outside the intended scope of this 
appendix. 
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SUB-APPENDIX B-l 


ERROR ANALYSIS OF STEADY STATE SOLUTION FOR A CONSTANT 
PROPERTIES SEMI-INFINITE RECEDING SOLID WITH CONSTANT 
SURFACE TEMPERATURE AND CONSTANT SURFACE RECESSION RATE 

B-l.l INTRODUCTION 

The differential equation for this problem is 

dt(x,9) d 2 t(x ) 0) ot(x 6) 

pc P — SS ’ * ~ £5 + P C P s — 5T— (B-i-i) 

Here the origin of the x-coordinate is fixed to the receding surface. Thus 
the left-hand side represents rate of energy storage (at constant x) • the 
first term on the right-hand side represents net heat conduction, and the 
final term represents net energy convection with the solid material. 

The associated boundary conditions are 

t (0, 0) - t e (B-l-2) 


t(x G) 


X — » oo 


« t 




It can readily be shown that a steady state solution of the system 
(B— 1-1) , (B-l-2), and B-l-3) is 



e - Sx/a 


(B-l-4) 


B-l. 2 FINITE DIFFERENCE FORM 

The Version 2 program uses Equation (B-l) as the finite difference form 
of the energy equation. in the steady state this relation becomes (consider- 
ing for convenience constant areas, constant properties, and constant nodal 
thicknesses) 


kA 9 
, 2 


”T' - 2T' + T* "| 

n-i n n+i 


+ pc S ^ 
h p d 


The difference forms in this equation have the 
expans ions : 


■ < T n + i " T A> " 0 

following associated Taylor 
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T (x + d) - 2 T (x) + T(x-d) 


+ * d T ( x ) + o (d 4 ) 

dx 2 12 ax 4 


E x ( B-l-6) 


T (x + d) - T (x) 
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Truncating the Taylor series expansions after the second term in both cases 
gives, for the two expansions, 


d T (x) _d" 

dx ? 12 dx 4 


(B-l- 8 ) 


where 


dT(x) . d d T( V 

dx 2 dx* 


d * x. ^ x + 


(B-l- 9 ) 


(B-l- 10 ) 


and 


- x_ £ x + 


(B— 1-11) 


Now the original differential equation for the steady state profile is 

( B- 1-12) 


kd t x d t 

— — + pc S — » 0 

h p dx 


dx 


or 


d 8 t + ( s.i dt 
dx ? dx 


( B-l- 13 ) 


From Equation (B-l- 13 ) may be derived the useful relation 

■ -in 


dj 

dx" 


dt 

dx 


( B-l- 14 ) 


Now if t m t, these relations may be employed ta simplify the expansions. 
Substitution of (B— 1 — 14 ) into ( B— 1 — 8 ) and (B— 1 — 13 ) into (B— 1 — 9 ) gives 


d 2 T(x) 

dx 2 


d^fsp dT(x i } 

12 l a I 


dx 


( B-l- 15 ) 


dT (x) _ d S\ dT (xg ) 

dx "" 2 iaJ dx 
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B-l . 3 ESTIMATION OF ERROR 


Equations (B-l-5), (B-l-7) , (B-l-7), (B-l-15) , and (B-l-16) imply that 

the solution to the difference equation corresponds to a solution of the 
differential relation 


kd 2 T 

5 + 

dx 2 


* dT 
pc S “ 
h p dx 



I ] 3 dT(x x ) 
a / dx 
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pc S 
K P 
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2 



dT (x ) 
2 

dx 


0 (B-l-17) 


If d is small, then « x, and so an approximate solution to this 

equation and its associated boundary conditions is 


T 



(B-l-18) 


where 
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Since 
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(B-l-21) 

so that the approximate 

solution is nearly 
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Since e~^ x / a i s the true solution, the relative error of the approximate 
solution is Bx/k, or 
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For the relative error to be small, the parameter Sd/* must be small. In 
this case, the relative error is given by 
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( B-l-24) 
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It is interesting to note that this expression very closely predicts the 
error magnitudes for the three cases shown in Figure B-8.* 

B-1.4 CONCLUSION 

The discussion has shown that errors in the steady .state profile of the 
constant properties semi-infinite receding solid with constant surface tem- 
perature and constant recession rate are given by Equation (B-l-24), provided 
that total errors, including those due to a poor choice of node sizes (con- 
duction errors), are small. 

The criterion for errors to be small is Sd/a « 1. This criterion 

pertains only to this problem, since it is a function of the exact solution. 
Since many ablation problems model this sample problem closely, the criterion 
should be of more general utility. 

In effect, Equation (B-l-24) gives an estimate of the size of errors due 
to the relatively crude convection term (errors 0(d)), provided that errors 
due to the more accurate conduction term (errors 0(d 2 )) are small. If con- 
duction errors are not small, errors due to the convection term are not 
readily estimated. Problems of estimating conduction term errors lie out- 
side the scope of this appendix. 


*Only in cases 2 and 3 was the nodal spacing sufficiently uniform to permit 
a true comparison. 
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APPENDIX C 


STUDY OF ALTERNATIVE TREATMENTS OF THE 
DECOMPOSITION KINETICS EQUATION IN 
THE CHARRING MATERIAL ABLATION PROGRAM 
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STUDY OF ALTERNATIVE TREATMENTS OF THE 
DECOMPOSITION KINETICS EQUATION IN 
THE CHARRING MATERIAL ABLATION PROGRAM 


C.l INTRODUCTION 

This appendix reports the results of a number of computational experi- 
ments aimed at defining and eliminating causes of irregular or uneven internal 
response solutions occasionally produced by the Aerotherm Charring Material 
Ablation Program (CMA) . For economy reasons, the calculations described here 
were limited to one material and a small number of transient conditions. 

Hence, it has not been possible to extract many general conclusions. However, 
the key problems have been fairly clearly identified, and in the process of 
doing so some significant improvements were made in the program operation. 

The appendix includes a number of graphs which illustrate the effect 
of various computational technigues* To present all the interesting aspects 
of the results would require an enormous mass of detail and many pages of com- 
puter output. Obviously a short memorandum can only touch upon a limited number 
of aspects, and hence the report necessarily omits many details and takes on 
a descriptive nature. 

C . 2 SUMMARY OF CONCLUSIONS 

The history of the instantaneous pyrolysis gas rate provides the most 
sensitive index to the smoothness of the internal response solution. Plots of 
this gas rate for slowly decomposing materials such as the phenolics are smooth 
and regular regardless of the transient. The solutions to certain rather rare 
problems have irregular gas rates. These do not affect the important overall 
aspects of the solution and so generally even these rare irregular solutions 
can be accepted. The boundary layer program (BLIMP) cannot accept irregular 
solutions, however. Efforts to eliminate oscillations for this application 
showed that: 

1. Replacement of the explicit treatment of density in the decomposition 
kinetic equation by exact integrated relations smoothed most solutions 
and shortened computation time 

2. Remaining unsmooth solutions could be smoothed by reductions in node 
size 

3. The number of nodelets per node has no importance for such problems. 
Nodelets may be eliminated, or greatly reduced in number, with great 
savings in computation time. 
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c. 3 general background 

Determination of the internal thermal response of an ablating, charring 
material requires the solution of a differential energy transport equation in 
depth. One form of this equation is 
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Here the origin of the x-coordinate is tied to the moving heated surface, 
whereas the origin of the y-coordinate is fixed in space (at the original loca- 
tion of the heated surface) . 

This energy equation is coupled to some decomposition kinetic equation 
for (3p/d0) . In the CMA program this decomposition equation is 


where 
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If all computations are done in terms of a nodal network tied to the 
receding surface, then to keep track of nodal (constant x) densities it is 
necessary to add to equations (C-l) and (C-2) above the convective relation 



Computational details for this equation will be discussed later. 

These three equations, with appropriate boundary and initial conditions, 
usually require finite-differencing for solution. Finite differencing may be 
done in a variety of ways, each of which has special aspects. 

The aspects of interest here pertain to the treatment of the decomposition 
equation (C-2) . If not treated properly in relation to the treatment of the 
energy equation (C— 1) , this equation can cause great computational difficul- 
ties. These difficulties are illustrated below. 
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C.4 FULLY EXPLICIT TREATMENT OF ENERGY AND DECOMPOSITION EQUATIONS, NODES ONLY 

Explicit" denotes that each "new" quantity, such as a nodal temperature, 
is computed in terms of "old" quantities only. Thus temperatures valid at the 
end of a time step are computed with temperatures and other properties values 
valid at the beginning of the time step. Similarly, a new density is computed 
from Equation (C-2) and (C-3) as 



(C-4) 


Here n refers to a fixed "node" or zone or lump of material defined by the 
finite differencing process. 

An important feature of explicit techniques is the instability of the 
finite-difference solution if the time step size exceeds some critical value. 

Even for time steps below this critical value, however, the solution for typical 
charring materials is apt to have severe oscillations due to irregularities in 
the decomposition calculation. This in turn is due to the extremely sensitive 
dependence of the decomposition rate, Equation (C-2), on both nodal density, since 
the reaction order m^ is usually at least unity, and in the nodal temperature, 
through the exponential term. For example, an increase (or an error) in tem- 
perature of 50° from one time step to the next can double or triple the decom- 
position rate. This results in excessive cooling of the nodes nearer the sur- 
face due to pyrolysis gas blow-by. These effects result in too low nodal 
temperatures and very much reduced decomposition rates. In the following time 
step, too much conduction occurs, resulting in too high temperatures again, and 
the cycle repeats. 

Irregularities of this type indicate that the decomposition event must be 
handled carefully. One such treatment is described in the next section. 


C.5 FULLY EXPLICIT TREATMENT OF ENERGY AND DECOMPOSITION EQUATIONS WITH NODES 
FOR THE ENERGY EQUATION AND NODELETS FOR THE DECOMPOSITION EQUATION 

One computational strategem to reduce this oscillatory decomposition 
problem is to employ a finer difference network for computing the decomposition 
than for computing the energy equation, with the aim of better characterizing 
the local decomposition temperatures and of more sharply defining the very 
steep density profile in the decomposition reaction zone. (One wishes to avoid 
this fine network for the difference form of the energy equation since the 
maximum time step allowed for stability reasons varies as the square of the 
nodal thickness.) 



This sub-network of fine nodes ("nodelets") has been used with good 
results. Experience shows that about ten nodelets per node adequately 
characterize the density profile and also damp out grossly oscillatory behavior 
for most materials of interest. 

Computation time with this method, although generally "acceptable", is 
long due to the large number of decomposition calculations required by the fine 
nodelet network. Nevertheless, an Aerotherm program based upon this procedure 
has proved to be of great utility in a wide variety of rocket nozzle and heat 
shielding calculations (Reference C-l) . The computation time of the program 
depends, of course, on the problem considered, but a typical problem requires 
three times real time (IBM 7094) . The next section describes a technique for 
shortening this time. 


C.6 FULLY IMPLICIT TREATMENT OF ENERGY EQUATION WITH NODES, FULLY EXPLICIT 
TREATMENT OF DECOMPOSITION AND MASS CONVECTION WITH NODELETS 

It has long been known that certain "implicit" treatments of the heat 
conduction equation, in which "new" temperatures are computed in terms of "new" 
temperatures, always give stable solutions regardless of the length of time step 
employed. For example, in simple transient conduction for constant material 
properties, a nodal equation would be 
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This formulation allows much longer time steps than the corresponding explicit 
expression . 

The more complex energy equation for a charring, ablating material, Equa- 
tion (C-l) , can be written in analogous fashion. The decomposition term op/o£) 
and the pyrolysis gas flow m , however, depend upon the nonlinear decomposition 
Equation (C— 2) , and render solution of the resulting series of equation very com- 
plex. The most straightforward solution method would be some iteration procedure, 
but it is not clear whether or not much time can be saved over the explicit 
procedure in this manner, since the number of iterations required may be large 
and various convergence testing techniques might have to be added. 

An alternative procedure is to leave both the temperature and density 
dependence decomposition events in explicit form, so that these nonlinear equa- 
tions can be solved explicitly. These results, in terms of dp/cW) y a ^d m g , 
can then be used in an implicit form of the iteration procedure; the suggested 
computational step implies that this one step gives adequate accuracy. 
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An earlier version of the present program incorporated this general 
scheme. This experience has shown that after the program worked "well" in the 
sense that the solutions usually displayed a regular and smooth behavior in 
contrast to the irregular or oscillatory behavior which one might expect for 
long time steps. For these successful problems, the time steps were much 
longer than those required in the analogous explicit solution and typical 
computation time was cut by a factor of three from that for the explicit pro- 
gram. Thus these problems ran in roughly real time (IBM 7094) . (Improvements 
to be described below will reduce this time still further.) 

Figure C-l shows a typical transient rocket nozzle solution. Although 
the example of Figure C-l indicates that this computational technique works 
well for most problems, experience has shown that the implicit treatment of 
the energy equation combined with a fully explicit treatment of the decomposi- 
tion equation does not always produce smooth solutions, particularly for rapid 
decomposers, such as nylon phenolic, exposed to large scale slow transients 
(a term to be defined later) , such as occur in reentry problems. 

Figure C-2 shows a sample calculation of such a problem in terms of the 
pyrolysis gas flow rate out the surface as a function of time. (This is the 
most sensitive indicator of the smoothness of the solution.) The material is 
a low density nylon phenolic. Decomposition reaction kinetics are computed 
with a three component model. Boundary conditions are a specified surface 
temperature beginning at 540° R and rising linearly in time to 5000° R at 10 sec- 
onds (and remaining at 5000 R thereafter) , combined with a specified surface 
recession of 30 mils/sec beginning at 10 seconds. Appendix A gives further 
problem specification details. 

Figure C-2 shows that during the ramp temperature input, before surface 
recession begins (at 10 seconds) , the computation has great difficulties in 
producing a smooth gas rate and displays both small irregularities and a large 
scale waviness. After surface recession begins the gas rate displays sharper 
oscillations of a smaller magnitude. 

To examine the effect of surface recession, the problem was rerun with 
c-he surface recession turn on at 5 seconds rather than 10 seconds. Figure c — A 
shows the pyrolysis gas rate results. 

The results display a number of very irregular appearing peaks and valleys 
The rough appearance of the solution is accentuated in the plot by the straight 
lines drawn between output points; since computation has in fact produced a 
number of intermediate points not shown, the actual course of the calculation 
may have been much smoother than Figure C-3 indicates. In any case, however, 
the gas rate is changing value in a rapid and presumably unrealistic manner. 



Fortunately, the uneven behavior shown in Figures C-2 and C-3 is not usually 
important. Abundant computational evidence, not to be reported here, indicates 
that the overall aspects of the solution (total surface recession, total pyroly- 
sis gas generation, surface temperature behavior, thermocouple temperature 
behavior, char thickness, and isotherm penetration) are affected in only a very 
minor degree by the irregularities. 

There is at least one important case, however, for which unevenness in 
the gas rate would be unacceptable, and that is for the coupled boundary 
layer solution (CABLE) . Convergence of the boundary layer computation will 
oe adversely affected by too large jumps in the pyrolysis gas injection rate. 

Therefore although uneven gas rate results are rare (being confined to a 
limited number of materials exposed to a special kind of transient) and not 
usually of importance, their influence on the boundary layer solution dictates 
a closer study attempting to eliminate uneven behavior. Such a study has been 
made and has resulted in significant improvements in computation procedures, as 
described below. 

C.7 TIME STEP LIMITATION BASED ON DENSITY CHANGE RATE 

An early device employed to smooth out irregularities limited the time 
step size according to an empirically determined limit on the density change 
per time step of the fastest decomposing component in the material. This device 
worked well in many problems, but for difficult cases it caused a sizable in- 
crease in computation time. Thus there was a strong economic incentive to re- 
move this "decomposition clamp". 

C 8 CHANGE IN TREATMENT OF DENSITY POTENTIAL IN DECOMPOSITION EQUATION FROM 
FULLY EXPLICIT FINITE DIFFERENCE FORM TO EXACT INTEGRATED FORM 

C.8.1 Introduction 

Up to this point the decomposition relation, Equation (C-2), and the con- 
vection relation, Equation (C-3) , have been treated in the following explicit 
manner: 

1. Given " old" temperatures and "old" nodelet desnities for each compo- 
nent i, "temporary" nodelet densities are computed with the finite 
difference convection relation (C-5) . 

2. These "temporary" densities are inserted into the decomposition 
kinetic relation (C-4) and "new" nodelet densities computed . 


C-6 



3. Pyrolysis gas evolution rates are computed from accumulated density 
change rates and are then used in the energy equation solution for 
" new" temperatures. 

This general procedure is a common idea and has been used previously (Reference 
(C-3) . 

Reference C-4 has pointed out, however, that the variables in the 
decomposition relation (C-2) , if temperature T is treated explicitly as 
is done here (T fixed during the integration), may be separated ana the 
integration on density performed exactly. The equation 
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can be written as 
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which for 7 1 integrates to 
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giving the density change rate 
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which yields the density change rate 
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This technique will obviously be "better" than the previous explicit 
treatment of the density, since a study of the computer results for the prob- 
lem shown in Figure C-3 indicates that much of the oscillatory behavior of the 
solution is due to the coupling between the decomposition and the energy equa- 
tions, which causes temperatures to be depressed after excessive decomposition, 
shutting off decomposition in the next step and allowing too great a tempera- 
ture rise again the the following step. This problem is now presumably elimin- 
ated and longer time steps should be allowable. In particular, the " decomposi- 
t ion clamp” should no longer be necessary. 

In fact, the improvement offer by the density integration procedure is 
very great, as will be seen in the examples cited below. 


C.8.2 Favorable Examples 

Figure C-3 provides graphic evidence of the powerful damping effect pro- 
vided by the exact integration forms of the decomposition kinetic equations. 

The figure displays the same problem as treated in Figure C-2 . Almost all the 
"high frequency" irregularities have been removed, and the computation steps 
much reduced (144 time steps at 11.9 seconds instead of 224 time steps in the 
older technique) . This saving is primarily due to the removal of the 
decomposition clamp. 

Figure C-5 shows the same problem as Figure C-4 but with the new exact 
decomposition relations. Here the smoothing effect is even more dramatic, 
especially in the region after the start of surface recession (5 seconds). 
Computation has been cut from 645 time steps to 171 time steps at 12 seconds. 

C.8.3 Less Favorable Examples 

Figures C-3 and C-5, when compared to Figures C-2 and C-4, indicate that 
a significant improvement has been made in computing technique. However, these 
same figures suggest that all is not yet well, and that the new decomposition 
relations are not a cure-all. In particular, in Figure C-3, the time between 
3 and 5 seconds (before surface recession starts) still shows a fairly unsmooth 
gas rate. Furthermore, there is a provocative change in solution character after 
10 seconds, when the surface temperature rise flattens out at a fixed value of 
5000°R. 3eth these aspects seem associated with tne same general pnenomenon: the 
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approach to steady state. Before 5 seconds there is no surface recession and the 
imposed linear surface temperature ri'e makes the problem everywhere "transient" , 
Between 5 and 10 seconds the imposed high surface recession rate greatly 
steepens the temperature profile and causes the problem to be more like a 
steady state problem. After 10 seconds the constant surface temperature allows 
the internal response to approach a steady condition very rapidly. 

Therefore, it appears that the computational technique is tending to have 
difficulties when the transient nature of the problem is dominating, that is, 
when the temperature profile and the pyrolysis zone (region of density change) 
are progressing through the nodal network; unless the decomposition is slow 
and the pyrolysis zone smeared out over many nodes (as in the case with the 
rocket nozzle problem shown in Figure C-l) . 

A slower and less artificial transient than that of Figures C-2 through C-5 
provides a good test of these suspicions. The problem represents the same 
nylon phenolic 3-component material, this time exposed to a ballistic entry. 

The material is exposed to air with a recovery enthalpy of about 15,000 Btu/lb; 
the heat transfer coefficient peaks at about 80 seconds. The first node is 
0.020 inches thick, 15 following nodes are 0.040 inches thick. Additional 
thicker nodes follow. (Complete specification of the problem conditions would, 
of course, be very lengthy. Since the problem serves only illustrative pur- 
poses here, further description is omitted (except for the property data in 
Appendix 3-1) .) Figure C-6 shows the pyrolysis gas rate as a function of time. 

The figure confirms the idea that a continuing transient condition can 
cause computational difficulties even for the improved decomposition procedure. 

An examination of the detailed nodal output reveals that in the period 0 to 70 
seconds the pyrolysis zone moves steadily through the nodal network from the 
surface node to approximately node 10. After 70 seconds the pyrolysis zone 
location stabilizes in the neighborhood of node 10. Figure C-6 shows that during 
the transient portion the surface pyrolysis gas rate is uneven (again it is per- 
tinent that the output plot does not show every point computed so that the actual 
succession of gas rate points would presumably have a smoother appearance) 
whereas after 70 seconds the output is much more smooth. 

A careful examination of other output information for this problem reveals 
that the unevennesses in the gas rate stem from a problem different from the 
coupled, on-off "oscillatory" problem observed before the introduction of the 
exact integrated decomposition relations. For example. Figure C-7 shows the den- 
sity histories of a number of nodes. This plot reveals two undesirable features 
of the computation: 

1. There is a tendency for one node to finish decomposing a little before 
the next node enters the rapid decomposition part of its history, so 
that for most of the time the surface gas rate is mainly determined 
by the events at a single node 



2, The progress of decomposition of each node, although apparently 
smooth, attains a maximum value when the nodal density is about 
midway between the virgin material density and the char density. 

Clearly these two aspects define the uneven appearance of Figure C-6, where the 
number of peaks corresponds closely to the number of nodes through which the 
pyrolysis zone has moved (this will later turn out not to be a general phenomenon, 
however) . 

Thus the unevenness problem seems associated with the nodal structure. 

This aspect is somewhat surprising, since the assignment of 10 nodelets per 
node was previously thought ample to define the pyrolysis zone and to provide, 
by linear interpolation between nodes, sufficiently accurate temperatures for 
decomposition. Indeed, a study of the node to node density variation (not to 
be presented here) indicates that the density profile is well defined by such 
a nodelet spacing for the node size used here. 

Additional computational tests were made to confirm the identification 
of the unevenness problem with the nodal structure. Figure C-8 shows a short 
calculation on essentially the same problem with the time step limited to 0.05 
seconds between 18 and 22 seconds and output called for every time step in the 
interval. (To save computation time the kinetic data has been changed to a 
single component model and the problem definition has been changed to one of 
specified surface temperature and recession rate. These changes cause no im- 
portant differences in the internal response.) The first plot on Figure C-8 
shows the gas rate plotted, as previously, every second. The second shows the 
gas rate between 18 and 22 seconds on an expanded scale, with every time step 
indicated. Study of the nodal density data (not presented here) shows that this 
expanded plot displays chiefly the decomposition of the second node. Figure C-8 
thus seems to confirm that the decomposition computation is smooth enough but 
tends to occur in waves associated primarily with the nodal size. 

This preliminary conclusion suggests a number of further diagnostic cal- 
culations aimed at understanding the exact nature of the unevenness in the gas 
rate and illuminating the role of nodal size. The sections below describe some 
of these calculations. 

C . 9 CHANGE FROM LINEAR INTERPOLATION BETWEEN NODAL TEMPERATURES TO CUBIC 
CURVE FIT 

To ascertain whether or not the linear interpolation between nodal tem- 
peratures to determine nodelet decomposition temperatures was perhaps imposing 
a nodal-scale inaccuracy on the total gas rates, the program was temporarily 
modified to employ a cubic curve fit between each pair of nodal temperatures 
to obtain more accurate temperatures for decomposition. The Aerotherm curve fit 
routines SL0PQ and 0GLE (Reference C-2) were used to generate the necessary 
parameters. Only minor improvements in gas rate smoothness resulted. 
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C.10 EXPERIMENTS ON NODAL SIZE 


C.10.1 Introduction 

The failure of a better interpolation scheme to produce smoother results 
points directly to nodal size as the dominant aspect, which in turn suggests 
that the uneven gas rate problem is directly connected to the accuracy of the 
energy equation solution. Additionally, this failure calls into question the 
utility of the nodelet concept in the present computation technique. Originally, 
nodelets were introduced to better define dcomposition events without having to 
accept any penalties in computation time in the old explicit solution. This 
computational technique proved to be extremely useful for ’'mild decomposers" 
such as the phenol ics, which decompose slowly over a range of temperatures. 

Little experience was gained with the explicit program on such materials as 
nylon,, which decompose very rapidly in a very narrow, sharply defined tempera- 
ture range, but even the brief experience indicated that computational perform- 
mance, in terms of a smooth gas evolution rate, was going to be poor. Now the 
more extensive experience described here with the implicit program confirms the 
difficulty of the "rapid decomposer" problem and suggests that the nodelet con- 
cept is less useful for this kind of problem. 

In all the calculations reported up to this point, there have been 10 
nodelets per node and the nodal sizes have been chosen with only one criterion 
in mind: to provide "enough" nodal points to define the temperature profile 

accurately. With a material like carbon phenolic and a typical transient, this 
procedure leads to temperature and density profiles roughly like: 



With nylon phenolics, for example, such a nodal sizing procedure is likely 
to yield temperatures and density profiles like: 
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Temperature 



Distance, with nodal 
centers indicated 


The use of nodelets in this problem thus serves to help define the sharp density 
profile, but the overall temperature accuracy is not likely to be sufficiently 
accurate to give meaningful instantaneous gas rates. 

This in fact appears to be the root of the problem of the uneven gas rate 
predictions in the examples considered earlier. The decomposition rate of nylon 
depends sensitively on the temperature and even small errors in nodal tempera- 
tures will cause wide excursions in the gas rate. For example, in the single 
component model for nylon phenolic, the decomposition equation is 
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where experimentally determined values of the constants are (References C-5 and C-6) 
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Therefore, when the decomposition rate is of importance it will be roughly 
proportional to 
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Thus a relative error in temperature of e causes a relative error in decom- 
position rate of e 30€ . For T = 1500°R, errors of 5°R, 10°R, 20°R, and 30° R 
give the following errors in decomposition rate 
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Hence even a 5° error in a nodal temperature will cause a 10% error in the gas 
rate, when the temperature is high enough for the decomposition to be important. 

A 5° error represents about 1/2% of the total nodal temperature rise in typical 
problems . 

The main drift of these observations then is that the accuracy of the 
temperature solution seems directly involved with the uneven gas rate problem. 

The following section adduces further evidence to this effect. 

C.10.2 Experiments with Smaller Node Sizes 

If accuracy of the temperature solution underlies the gas rate problem, 
the only solution to the trouble consists of smaller node sizes. Figure C-9 
shows the same problem as illustrated earlier in Figure C-6 (aside from the minor 
changes to a single component decomposition and a specified surface temperature 
and recession rate problem rather than a surface thermochemistry problem) , 
but with 0.02 inch nodes instead of the former 0.04 inch nodes. The improve- 
ment is striking, the amplitude of the gas rate "waves” being reduced by about 
half. 

Thus it seems that nodal size is the dominant parameter, and that perhaps 
the number of nodelets per node is more or less irrelevant for this particular 
problem. To check this aspect, the problem of Figure C-9 was recomputed once 
again with 0.008 inch nodes, but with only two nodelets per node. This yields 
the same nodelet size as for the problem of Figure C-6. The results, shown in 
Figure C-10, may be compared to the earlier results for the same problem shown 
in Figures C-6 and C-9. Again the "wave amplitude” has been much reduced. 
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Figure C-10 has a conspicuous dip in the gas rate at 25 seconds. It 
happens that the specified surface temperature boundary condition for this 
problem has a marked decrease in slope at 25 seconds; to check if perhaps this 
caused the gas rate dip the problem was rerun with a ramp input for surface 
temperature and recession rate. Figure C- 11 shows the results and proves that 
irregular behavior is characteristic of the computational system rather than 
the specific boundary condition. 

As a final test of the effect of smaller nodes on the smoothness of the 
solution, the ramp input problem of Figure C-ll was rerun with a new nodal 
size distribution. The first node was 0.008 inch, and the next seven nodes 
were 0.016 inches thick. These nodes were followed by a- number of nodes of 
0.004 inches thickness, half as large as the nodes of the problems illustrated 
in Figures C-10 and C-ll. Figure C-12 shows the results of this last computa- 
tion. As expected, during the time the pyrolysis zone is progressing through 
the thick nodes (in this problem 0 to 30 seconds) the gas rate is very 
irregular. As soon as the pyrolysis zone passes into the thin nodes, the gas 
rate smooths out almost completely, except for one small wave at 38 seconds. 

The computations illustrated in Figures C-6, C-9, C-10, C-ll, and C-12 
thus indicate that the choice of nodal size and number of nodelets per node 
determines the smoothness of the pyrolysis gas rate. It would, of course, be 
desirable to have quantitative criteria for selecting these quantities, but 
the overall complexity of this problem seems to dictate a trial and error 
procedure. To facilitate the selection of these quantities, Appendix C-2 pro- 
vides some useful estimates of program execution time as a function of the 
number of nodes and nodelets. 

C.ll CONCLUSIONS 

The early versions of the Aerotherm Charring Material Ablation program, 
which treated the decomposition kinetic equations explicitly in both density 
and temperature and which used an implicit formula for the energy equation in 
depth, produce excellent results for slow decomposers such as phenolics 
regardless of the type of transient involved. Fast decomposers such as nylon 
cause irregularities in the solution values of quantities of interest, 
especially in the rate of pyrolysis gas generation, for those particular 
transient problems in which the pyrolysis zone is passing rapidly through the 
nodal network (which is tied to the receding heated surface) . 

Usually these irregularities cause no problems since experience shows 
that they do not affect the important overall aspects of the solution, such 
as total surface recession and char thickness, for example. However, since 
it is expected that irregularities in the pyrolysis gas rate will cause 
convergence difficulties in the coupled boundary layer program (BLIMP) , the 
present study included attempts to eliminate irregular or wavy gas rate 
predictions. 
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An early computational strategem lumped the rapid decomposition effect 
into one of the available three decomposing constituents and then applied a 
time step limit to the overall solution to prevent the decomposition rate 
(density change rate) of this component from exceeding some empirically deter- 
mined amount. This semi-rational time step limit was quite effective, but 
turned out to be rather costly for nylon, which for many transients requires very 
short time steps for a smooth solution. 

A better smoothing procedure abandons the explicit treatment of density 
in the decomposition relations in favor of exact integrated relations for the 
density dependence. The implicit-density nature of these relations exerts a 
powerful smoothing effect on the solution, and allows removal of the time step 
clamp on decomposition rates with a consequent saving in the number of solution 
steps. Although the integrated decomposition relations are about twice as slow 
to compute as the simpler explicit relations, the saving in time steps in diffi- 
cult to handle problems was observed to be a factor of four to five in the num- 
ber of time steps. Since the computation time is mainly devoted to decomposi- 
tion, net computation time has been cut by about a factor of two. 

Despite the strong smoothing effect of the implicit density relations, 
certain particularly unfavorable combinations of materials and boundary conditions 
can still cause markedly uneven gas rate predictions. The roots of this behavior 
could not be definitely identified due to the complexity of the internal response 
solution. Computational experiments demonstrated, however, that nodal size is 
the key parameter affecting the smoothness of the solution, which points to the 
overall accuracy of the energy solution as the key aspect. Smaller node sizes 
resulted in much smoother solutions and no doubt remains that entirely smooth 
solutions can be obtained with "sufficiently" small nodes. With small nodes 
the nodelet structure becomes of minor importance, and the number of nodelets 
can be cut down to the minimum (within the present program logic) of two without 
harm to the solution. Thus, the net computation time is in general not increased 
by the smaller node size. 

It would, of course, be desirable to have a quantitative criterion for 
selecting nodal and nodelet sizes, but the overall complexity of the problem 
appears to preclude a successful analysis and to demand an empirical approach. 

To date, computational experience has not been sufficient to define any general 
rules giving good nodal sizing. 
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■SUB-APPENDIX C-l 
NOTES ON PROBLEM SPECIFICATION 


C-l.l PROBLEM OF FIGURES C-2 TO C-5 

The material is a low density nylon phenolic with a resin mass fraction 
of 0.5, resin residual mass fraction of 0.5, and a nylon residual mass fraction 

of zero. Decomposition kinetic data, adapted from the results of Reference 9, 
were 


Component 

H o, 

l 

p r. 

1 

k. 

l 

Ea ./R 

l 


(lb/ft 3 ) 

(lb/ft 3 ) 

( sec” 1 ) 

(°R) 

A 

6.028 

0 

1.4 x 10 4 

15,400 

B 

18.084 

12.056 

4.48 x 10 9 

36,400 

C 

71.0 

0 

1.85 x 10 13 

47,700 


m . 

l 


3 

3 

1.5 


Enthalplies of formation, reference temperature 536°R: 

Virgin Plastic 0 

Char 0 

Pyrolysis gas - 1311 Btu/lb 

Material Properties: 


Virgin Plastic 


Char 


T(°R) 


Specific Heat 
(Btu/lb) 


Thermal 
Conductivity 
(Btu/ft sec ° F) 


460 

. 360 

1.5 x 10” E 

550 

.427 

1.39 

650 

,475 

1.39 

7 60 

.500 

1.39 

860 

.500 

1.62 

910 

.500 

1,68 

960 

.500 

1.69 

1560 

. 500 

1.69 

6460 

. 500 

1.69 

460 

. 122 

, 384 x 10”' 

1000 

. 305 

.650 

1500 

.405 

.900 

1960 

.480 

1.183 

2110 

. 500 

1.275 

2260 

.520 

1. 375 
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T ( J R) 


Specific Heat 
( Btu/lb) 


3260 
3360 
3660 
3760 
3860 
4060 
4160 
4310 
4610 
4760 
5460 
6460 

Pyrolysis gas sensible enthalpy: 

T(°R) 

900 
1800 
2700 
3600 
4500 
5400 

Nodal sizes: 

1 at 0.015 inches, 9 at 0. 

1 at 0.090 inches 

10 nodelets/node 

Boundary conditions are described 


Thermal 
Conductivity 
(Btu/ft sec °F ) 


. 600 

1.891 

.601 

1.947 

.605 

2. 109 

.606 

2. 156 

. 608 

2.212 

. 611 

2. 320 

.612 

2.416 

. 615 

2.539 

. 619 

2 . 966 

.621 

3.294 

.630 

4.931 

.630 

4.931 


h (Btu/lb) 

-1255 . 
10.0 
2687. 
3694. 
6010. 
12006. 


inches, 2 at 0.050 inches, 


in Section 6. 


C-1.2 PROBLEMS OF FIGURES C-6 TO C-12 

The problems had a low density nylon phenolic slightly different from 
that above. For the 3-component model, the kinetic coefficients were the same 
as those above except that m^ * 1. Densities were 
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Component 


A 

B 

C 


(lb/ft 3 ) 

6.55 

19.60 

71.00 


r . 

iwi tn 

0 

18.68 

0 


For the single component model the kinetic data is given in Section 10.1 above. 
Densities for this case were 


Component 

A 

B 

C 


(Ib/ff 

35.0 

0 

35.0 


r i 

db/ft 3 ) 

15.0 
0 

35.0 


Enthalpies of formation, reference temperature 536°R 


Virgin Plastic 
Char 

Pyrolysis gas 


0 

0 

-1117 Btu/lb 


Material Properties: 



T C R) 

Specific Heat 
(Btu/lb) 

Thermal 
Conduct ivit 1 
(Btu/ft sec ! 


460 

. 29 

1.26 x 10' 


660 

.44 

1.34 

Virgin Plastic 

760 

.50 

1.37 


860 

.54 

1.40 


960 

.55 

1.42 


1260 

.55 

1.55 


6000 

.55 

1.55 

Char 

500 

. 10 

.10 x 10’ 


1460 

. 39 

1.80 


1960 

.49 

2.46 


2460 

.50 

2.90 


3460 

.50 

4.00 


4460 

.50 

5.60 


4760 

.50 

6.40 


5500 

.50 

7.50 
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elastic omissivity was 0.9; char emissivity was 0.6. 


Pyrolysis gas sensible enthalpy: 


T(°R) 

900 

1800 

2500 

2700 

3000 

3600 

4500 

5400 


h^ (Btu/lb) 

-2206 
- 854 
459 
1021 
1967 
2757 
3841 
5850 


Boundary conditions: 


Ramp Transient 


Time ( sec) 

Lt 

T (°R) 
w 

S (mils/ 

C-10) 0 

540 

0 

5 

1430 

0 

10 

1590 

0.08 

17 

1820 

0.15 

22 

2170 

0.30 

25 

2350 

0.50 

30 

2850 

0.85 

35 

3600 

1.30 

40 

4100 

1.80 

45 

4500 

2.28 

50 

4770 

2. 75 

55 

5100 

3.30 

60 

5400 

4.25 

65 

5700 

5.80 

70 

5875 

7.50 

75 

5850 

8.00 

80 

5750 

8.42 

) o 

540 

0 

40 

4000 

1.0 
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SUB-APPENDIX C-2 


ESTIMATION OF EXECUTION TIME 


The following table of the number of operations in each section of the 
charring material ablation program provides a useful estimate of execution 
time as a function of the number of nodes, the number of nodelets per node, 
and so on. 

The input subroutine and all other operations, such as output, not per- 
formed each time step, have not been included. Only floating arithmetic is 
considered . 


Table of Floating Operations 
for One Time Step 



Add 

Subtract 

Multiply 

Divide 

Exp f 

Log f 

Mi seel laneous 
Nodal Computa- 
tions, per Node 

32 

36 

64 

25 

0 

0 

General Decompo- 
sition, per 
decomposing Node 

22 

15 

13 

3 

0 

0 

Component Decompo- 
stion, per Com- 
ponent for each 
Decomposing 
Nodelet* 

1 

3 

4 

3 

3 

2 

Surface Calcula- 
tions, per 
Iteration 

25 

39 

50 

in 

1 

0 


For the exact integrated relations and m ? 1. For the simpler explicit 
decomposition relations the numbers in this row were 1, 4, 4, 2, 1, 1. 


With 
multiply = 
yields the 


the average operating times for the IBM 7094 (add and subtract = 14 p 
7 ps } divide - 12 ps, exp f ** 188 ps, log f » 226 ^s) , this table 
following floating arithmetic time estimate per time step. 


£(1136 • I + 645) J . F + 1700 N + K(1534) M sec/step 



where 


T = execution time for a single time step (psec) 

I = number of decomposing components in the material 
j = number of nodelets per node 
F = fraction of nodes decomposing 
N = total number of nodes 

K - number of iterations in surface energy balance 

To estimate total execution time, T may be multiplied by the total num- 
ber ol Lime steps, a number approximately determined by the user in the choice 
uf the maximum time step allowed. F should be an average value for the frac- 
tion of nodes actively decomposing during the probLem {usually between 1/3 and 
1/2) . A study of numerous problems reveals that K averages very close to 3. 

The resulting total time must be multiplied by a factor to account for 
" a- Imi n i s t r. 1 1 ivc M arithmetic such as DO loop indexing. Experience with this 
program on the 7094 suggests a factor between 1.4 and 1. >. 

Additional machine time will be consumed by input and loading operations, 
overlay or special tape handling (if necessary), and output operations, the 
details of which are peculiar to individual computing facilities. 
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TRANSFORMATION OF THE IN-DEPTH ENERGY EQUATION 









APPENDIX D 

TRANSFORMATION OF THE IN-DEPTH ENERGY EQUATION 

The in-depth energy equation, which is Equation (38) in the main text 
above, is 


PC 


dT 


p at 


1 = T 

d 


- h 

Jx A 

dx | 

1 He 

p 6 . 


+ {ph) 9 s 


^2 Sh q | 

A Sx Je 


+ h 


* 2 .} 

g ae] 


y 


where 


h 


P p^p p c^c 
p p ~ 


(D— 1) 


(D-2) 


This form of the energy equation turns out to be convenient for differenc- 
ing and for machine treatment, but it lacks a certain amount of physical 
clarity, primarily because no single "energy of decomposition" term can be 
identified. Equation (D-l) can be cast into more appealing form by definitely 
identifying this quantity. The necessary operations involve chiefly an expan- 
sion of the a/Sx(ph) term in Equation (D-l), as follows; 

Since Equation (21) of the main text gives 


ph = € p h + (1 - e ) p h 

p K p p v p c c 


then, dropping the "constant-0" subscript for convenience temporarily. 


3LP.M. 

3x 


|- (e p h ) - |- (e p h ) + |- (p h ) 
dx p M p p' dx pc c' dx VH c c 


(0-3) 


ft ft 

(p h - p h ) -r — ^ + e 3— (oh - p h ) + -r— (p h ) 
PP cc dx p >dx PP *c c dx M cc 

(D-4) 


^ ^ ft m ft rn 

(p h - p h ) + e p C (1 - e ) P C - 7 ^ 

P P c c dx P P p dx v p h c p dx 

P c 


(D-5) 
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But by the defining Equation (28) of the main text 


pC = e p C + (1 - e ) p C 
M P P^P P p P p c 


(D-6) 


so that 


S(ph)_ 

dx 


dG 


(p h - p h ) - — * 
v ^p p K c cr dx 


, ^ dT 

+ pC T 

K P 


(D-7) 


Now, however, we have by Equation (20) of the main text that 


p = ep + ( 1 - e ) p 
p P p c 


(D-8) 


so that 


le. 

dx 


dc^ 


de_ 


p p dx P c dx 


(D-9) 


dx 


ip . 1 

dx , 
It 


P P " p c 


( D- 10) 


Substituting this into Equation (D-7) gives finally 

Ae.1 


1 

dx 


-Jo 


n Sx e p p dx 


(D-ll) 


where h is defined by Equation (D~2) above. 

Now substitution of Equation (D-ll) into the original energy Equation (D-l) 

gives 


_ dT 

pC p 9 0 


1 

d_ 

kA |2 

+ h 

s |fi.| 

dp 

^ X 

A 

X 

dx ! 

1 9x J 

9x| 0 

d 0 J 



+ SpC 

p a*j 0 


m dh ^ 

+ -2 

A dx 


+ h 


lfi.1 

g 96 j 


(D-12) 


D-2 



Hut the conservation of mass Equation (15) of the main text gives 


la 

3 G 



+ S 

5x 


0 


(D-13) 


Putting this into the h coefficient above and collecting terms we obtain 
the desired final relation 


pC 


3T 

3 6 


1 

A 


1_ 

3x 


kA 


3T I 
3x i 


(h g 


h) 


ll> | 

30 


/y 


A _ 3T 
SfjC 3 — 
p 3x 


m 

A 


oh 

3x 
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This equation may be viewed as a "conservation of sensible energy" equa- 
tion at constant x, and thus has a close relation with the usual form of 
energy equations in applications including chemical changes. The second term 
on the right hand side represents the "creation of sensible energy" due to 
pyrolysis . 
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